A simple coupled potential QM/MM/MD simulation.

(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

SECTION 4

A Coupled Potential QM/MM Simulation

By Ross Walker

Updated for AMBER 15

4) Comparing the results

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) =   455.28  PRESS =     0.0
 Etot   =     -9086.7258  EKtot   =      4144.0991  EPtot      =    -13230.8249
 BOND   =         0.0650  ANGLE   =         0.1616  DIHED      =         2.5068
 1-4 NB =         1.1157  1-4 EEL =       -19.4429  VDWAALS    =       930.8446
 EELEC  =    -14146.0758  EHBOND  =         0.0000  RESTRAINT  =         0.0000
 Ewald error estimate:   0.1330E-03
 ------------------------------------------------------------------------------

 NSTEP =        1   TIME(PS) =       0.002  TEMP(K) =   346.22  PRESS =     0.0
 Etot   =    -10079.4435  EKtot   =      3151.3814  EPtot      =    -13230.8249
 BOND   =         0.0650  ANGLE   =         0.1616  DIHED      =         2.5068
 1-4 NB =         1.1157  1-4 EEL =       -19.4429  VDWAALS    =       930.8446
 EELEC  =    -14146.0758  EHBOND  =         0.0000  RESTRAINT  =         0.0000
 Ewald error estimate:   0.1330E-03
 ------------------------------------------------------------------------------

QMMM

 NSTEP =        0   TIME(PS) =       0.000  TEMP(K) =   455.28  PRESS =     0.0
 Etot   =     -9070.8407  EKtot   =      4144.0991  EPtot      =    -13214.9397
 BOND   =         0.0000  ANGLE   =         0.0000  DIHED      =         0.0000
 1-4 NB =         0.0000  1-4 EEL =         0.0000  VDWAALS    =       901.0150
 EELEC  =    -14062.6280  EHBOND  =         0.0000  RESTRAINT  =         0.0000
 PM3ESCF=       -53.3267
 Ewald error estimate:   0.5644E-01
 ------------------------------------------------------------------------------
 NSTEP =        1   TIME(PS) =       0.002  TEMP(K) =   346.23  PRESS =     0.0
 Etot   =    -10063.4652  EKtot   =      3151.4746  EPtot      =    -13214.9397
 BOND   =         0.0000  ANGLE   =         0.0000  DIHED      =         0.0000
 1-4 NB =         0.0000  1-4 EEL =         0.0000  VDWAALS    =       901.0150
 EELEC  =    -14062.6280  EHBOND  =         0.0000  RESTRAINT  =         0.0000
 PM3ESCF=       -53.3267
 Ewald error estimate:   0.5644E-01
 ------------------------------------------------------------------------------

The first thing you should notice is that the QM/MM result has an extra field (PM3ESCF). 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 QM/MM 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 QM/MM mdcrd file into vmd:

NMA topology file (AMBER7 Parm ) and then load the md_qmmm.mdcrd file into this molecule. Remember, we have used periodic boundaries here so select "Amber Coordinates with Periodic Box".

We should again remove the water molecules so our NMA is easier to see.

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 QM/MM loaded trajectory to the beginning and then load a new molecule. Select the NMA.prmtop again and then load the 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 color the two molecules differently in Graphics->Representations if you want).

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


Classical (MM)	QM

I leave it to you the reader to try other forms of analysis, you could look at the difference between average structures, plot the O-C-N-H out of plane angle as a function of time. You could also try cluster analysis and hydrogen bond analysis. How does the frequency of hydrogen bonds differ between the classical and the QM/MM simulation? See tutorial B3 for details on this type of advanced analysis.

You could also try coupling this QM/MM approach with Nudged Elastic Band. A classical (MM) introduction to this is given in tutorial A5. By utilizing QM/MM it becomes possible to look at reaction pathways instead of just conformational changes.

Return to index