(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.)

TUTORIAL A2

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.

(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.)