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

TUTORIAL 6 - AMBER 9 VERSION - SECTION 2

A Coupled Potential QM/MM Simulation

By Ross Walker

PLEASE NOTE THIS TUTORIAL IS FOR AMBER 9 ONLY.

2) Classical MD

Now we have our topology and coordinate files we are ready to carry out classical molecular dynamics. First things first though we need to minimise our system to remove any bad contacts created by solvation and also to allow our artificially created NMA structure to relax. For this we will run 500 steps of minimisation. Our input file is very similar to files we have used in the previous tutorials. Since AMBER 9 now supports periodic boundaries with PME for both classical and QM/MM MD we will use it (NTB=1). Amber 9 also supports shake for QM atoms and since we are not interested in any kind of reaction involving hydrogen in the QM/MM simulation we will use shake (NTC=2, NTF=2) for this simulation to allow us to use a 2 fs time step. Since this is a periodic simulation with PME we can safely use a cut off of 8 angstroms. Compare this with AMBER 8 where we had to use a solvent cap and thus a cut off of 20 angstroms. Here is our input file:

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

Input Files: min_classical.in, NMA.prmtop, NMA.inpcrd
Output Files: min_classical.out, min_classical.rst

Once we have minimised the system we will run 1,000fs of molecular dynamics at a constant temperature of 300K. Normally we would slowly heat our system to 300K but since this is just a small NMA molecule in solution it should be fairly stable and should survive being started out at 300K. During the MD we will write to both our output file and mdcrd file every step. In this way we will be able to see the detailed motion of our system. At the expense of slowing the simulation down while it writes a considerable amount of information to disk. Here is the input file for the MD run:

Now we run the md.

$AMBERHOME/exe/sander -O -i md_classical.in -o md_classical.out -p NMA.prmtop -c min_classical.rst -r md_classical.rst -x md_classical.mdcrd

This should take about 3 mins to run.
Input Files: md_classical.in, NMA.prmtop, min_classical.rst
Output Files: md_classical.out, md_classical.rst, md_classical.mdcrd

Now, let's do two things. First of all we will visualise our MD trajectory in VMD and secondly we will use ptraj to create an average structure for us.

Step 1 - VMD. Fire up VMD and load the NMA topology file (remember to select parm7 [vmd1.8.3] or AMBER7 Parm [vmd1.8.4]) and then load the md_classical.mdcrd (remember that this coordinate file has box information in it) file into this molecule. I gzipped it to save space so make sure you decompress it first.

While it is good to see our big box of water "floating" about we are more interested in the dynamics of the NMA residue so let's remove the water.

Graphics->Representations
In Selected Atoms type:
all not water
And hit apply
At the same time change the drawing method to CPK.

You should now be able to watch a movie of our NMA molecule "wiggling" about.

The thing to notice is how planar the amide C-N-H group remains during the simulation. This is largely artificial since the nitrogen atom has a lone pair which acts to force the hydrogen out of plane. This is a deficiency of the classical force field as it is currently implemented. Work is under way to try to improve this.

Lets repeat the simulation with QM/MM.