Amber Workshop - tutorial 6 - AMBER 8 VERSION

(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

PLEASE NOTE THIS TUTORIAL IS FOR AMBER 8 ONLY.
AN AMBER 9 QM/MM TUTORIAL IS AVAILABLE HERE

We will now qualitatively compare the results of the two simulations. First of all lets take a brief look at the two output files:

Classical

 NSTEP =        0   TIME(PS) =       0.000  TEMP(K) =   305.74  PRESS =     0.0
 Etot   =     -1705.9816  EKtot   =      1127.3405  EPtot      =     -2833.3221
 BOND   =       344.6722  ANGLE   =         0.6147  DIHED      =         2.4711
 1-4 NB =         0.8346  1-4 EEL =       -19.9098  VDWAALS    =       935.3651
 EELEC  =     -4097.3700  EHBOND  =         0.0000  RESTRAINT  =         0.0000
 ------------------------------------------------------------------------------


 NSTEP =        1   TIME(PS) =       0.001  TEMP(K) =   305.75  PRESS =     0.0
 Etot   =     -1705.9601  EKtot   =      1127.3620  EPtot      =     -2833.3221
 BOND   =       344.6722  ANGLE   =         0.6147  DIHED      =         2.4711
 1-4 NB =         0.8346  1-4 EEL =       -19.9098  VDWAALS    =       935.3651
 EELEC  =     -4097.3700  EHBOND  =         0.0000  RESTRAINT  =         0.0000
 ------------------------------------------------------------------------------

QMMM

 NSTEP =        0   TIME(PS) =       0.000  TEMP(K) =   305.74  PRESS =     0.0
 Etot   =     -3525.1146  EKtot   =      1127.3405  EPtot      =     -4652.4551
 BOND   =       343.7005  ANGLE   =         0.0000  DIHED      =         0.0000
 1-4 NB =         0.0000  1-4 EEL =         0.0000  VDWAALS    =       982.1865
 EELEC  =     -5917.1502  EHBOND  =         0.0000  RESTRAINT  =         0.0000
 ESCF   =       -61.1920
 ------------------------------------------------------------------------------


 NSTEP =        1   TIME(PS) =       0.001  TEMP(K) =   305.75  PRESS =     0.0
 Etot   =     -3525.0931  EKtot   =      1127.3620  EPtot      =     -4652.4551
 BOND   =       343.7005  ANGLE   =         0.0000  DIHED      =         0.0000
 1-4 NB =         0.0000  1-4 EEL =         0.0000  VDWAALS    =       982.1865
 EELEC  =     -5917.1502  EHBOND  =         0.0000  RESTRAINT  =         0.0000
 ESCF   =       -61.1920
 ------------------------------------------------------------------------------

The first thing you should notice is that the QMMM result has an extra field (ESCF). This is the energy due to the quantum part of the calculation (the NMA) within the presence of the charge field of the MM part (the water). You should also notice here that the QMMM result lacks ANGLE, DIHEDRAL, 1-4 NB and 1-4 EEL energies. This is exactly as we expect since the NMA bonds, angles, dihedrals and VDW and electrostatic terms are now dealt with inside the QM calculation. The only atoms remaining in the MM section are TIP3P water molecules. These are actually triangulated water and so only have bond terms (TIP3P water does not have an angle component). The water molecules are also only 3 atoms in size and so do not have any 1-4 NB or 1-4 EEL interactions.

If you compare the two output files further you should find that the temperatures are reasonably stable and similar. They will not be exactly the same since we are tracking different trajectories.

Now, lets compare the actual dynamics. First off lets load the QMMM mdcrd file into vmd:

Fire up VMD and load the NMA topology file (remember to select parm7) and then load the NMA_md_qmmm.mdcrd file into this molecule. Remember, we have not used periodic boundaries here so select "crd" (not "crdbox").

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. How does the dynamics of the NMA differ from the classical simulation? If you want you can load both molecules into VMD and animate them together.

Rewind the current QMMM loaded trajectory to the beginning and then load a new molecule. Select the NMA.prmtop again and then load the NMA_md_classical.mdcrd file into this molecule. When you close the load molecule dialog box you should find that both molecules are displayed. You can now play through the trajectory slowly and watch the difference. (You can colour the two molecules differently in Graphics->Representations if you want).

You should notice that the dynamics of the NH group are different in the QMMM calculation. The nitrogen regularly forms a pyramidal structure and the amine hydrogen is more often out of the O-C-N plane. The QMMM behaviour here is believed to be the more accurate. The NH backbone dynamics are known to be troublesome to model accurately with classical force fields.