QM/MM tutorial

QM/MM calculations on a Diels-Alder antibody catalyst.

Introduction

We are now going to see how reaction barrier changes when we solve the reactants in water. Since water is a polar solvent we can anticipate that there will be an effect of some sort. In this part we first create and equilibrate a simulation box with the reactants and water molecules. Then we perform a Linear Transit calculation, using the reaction coordinate of the previous step (figure 3).

QM/MM subdivision

The system consists of the two reactant molecules solvated in water and one Na+ ion. The ion compensates the overall charge of -1 on the aliphatc tail (-R, figure 1)). This system is way too big to be treated at the QM level. Therefore we divide the system in a small QM part and a much bigger MM part. The QM part consists of the reactants, without the aliphatic tail and is described again at the semi-empirical PM3 level, while the remainder is modelled with the GROMOS96 forcefield. Figure 6 shows the subdivision used in this part of the tutorial.

In this part of the tutorial we are going to perform again a Linear Transit calculation, but this time, the reactants are fully solvated. The details on how to perform a QM/MM Linear transit calculations in gromacs were discussed previously.

back to top

Finding product, reactant and transition state geometries in water, using Linear Transit

We will again make use of scripts to perform the Linear Transit.

• create_tops.scr
• run.scr
• get_ener.scr

• We start again by creating a topology. In this case we only need to add to the topol.top a line stating how many waters there are, how many and what type of ion we have. We also need to add include statements for the water model and ion model we use:

  #include ions.itp
  #include flex_spc.itp

  A complete topology file should look like this topol.top, which is available for download.

• We then need an index file and an mdp file. The only difference with the files we used in the previous part of this tutorial is that now the water molecules and ion are included. We will not create them ourselves now, but use the ones available (index.ndx, LT.mdp)instead.

• Now that we have the input files and scripts, we will perform the Linear Transit calculation. This time, we let the reaction coordinate vary in two stages. First, we go from 0.15 to 0.4 nm in 200 steps. Then we do an additional calculation in which we go down from 0.15 to 0.12 in 50 steps. When both stages are done we put the outcomes together in one graph. The reason for performing two seperate Linear Transits, instead of only one, is that it is very hard to get a converged structure starting at 0.12 nm on the reaction coordinate. Much better convergence is achieved when starting at 0.15 nm, the value of the reaction coordinate in the starting structure(confin.gro) and slowly go down to 0.12 nm. Also note that we use a more loose convergence criterium for the minimization of the points (emtol = 600) in the LT.mdp file). The system is much bigger than before, and therefore harder too mimimize.
  • create a new subdirectory for the Linear Transit in water:

    localhost:~>mkdir LTwater

  • download the scripts, confin.gro, topology, index and LT.mdp into that directory

  • create two additional subdirectories in the LTwater water subdirectory, one for increasing and one for the decreasing the reaction coordinate:

    localhost:~>mkdir up

    localhost:~>mkdir down

  • go into the up subdir and execute the create_tops.scr scripts to create 201 subdirectories (called step_0, ..., step_200) and create a topol_A.itp with different constraints lengths in the range 0.15 to 0.4 nm.

    localhost:~>cd up

    localhost:~>../create_tops.scr 0.15 0.4 200

  • execute the run.scr scripts to run the minimizations in the different subdirectories (called step_0, ..., step_200). To speed up the convergence, the script takes the output coordinates of the previous Linear Transit point as input in the current minimization

    localhost:~>../run.scr 0.15 0.4 200

    This will take a while, depending on the speed of your computer system.

  • execute the get_ener.scr script to retrieve the energies from the individual linear transit points. Redirect the output to a file, let's call it eqmmm.xvg.

    localhost:~>../get_ener.scr 0.15 0.4 200 > eqmmm.xvg

  • Next we go back one subdir back and enter the down subdir to do the other linear transit.

    localhost:~>cd ../down

    execute the create_tops.scr scripts again to create 51 subdirectories (called step_0, ..., step_50) and create a topol_A.itp with different constraints lengths in the range 0.15 to 0.12 nm.

    localhost:~>../create_tops.scr 0.15 0.12 50

  • execute the run.scr scripts to run the minimizations in the different subdirectories (called step_0, ..., step_50). To speed up the convergence, the script takes the output coordinates of the previous Linear Transit point as input in the current minimization

    localhost:~>../run.scr 0.15 0.12 50

    This will also take a while.

  • execute the get_ener.scr script to retrieve the energies from the individual linear transit points. Redirect the output to a file, let's call it eqmmm.xvg.

    localhost:~>../get_ener.scr 0.15 0.12 50 > eqmmm.xvg

• We now have a look at the result of the Linear Transit.
  • Go back one subdir and import up/eqmmm.xvg and down/eqmmm.xvg in xmgr to see the energy as a function of the reaction coordinate.

    localhost:~>xmgr down/eqmmm.xvg up/eqmmm.xvg

    

    The maximum is at 0.20875 nm (-1062.9 kJ/mol). This corresponds to point 47 in the up subdirectory. Let's have a look at that structure. We simply go into the directory up/step_47 and use editconf to convert the confout.gro into confout.pdb

    localhost:~>editconf -f confout.gro -o confout.pdb

    We then use rasmol to visualize the structure:

    localhost:~>rasmol confout.pdb

    

    We show the first shell of water.

    • RasMol>restrict within(6.0,1)
    • RasMol>zoom 600
    • RasMol>wireframe 40

    Clearly there are many water molecules interacting with the molecule.

    The minima are at 0.335 nm (-1202.4 kJ/mol) and 0.15125 nm (-1271.7 kJ/mol) and correspond to the reactant and product states respectively.

    Now, we plot the energy curve we obtained previously for the reaction in vacuo in the same figure to see the effect of the water on the energetics of the reaction. We modify a bit the offset of the vacuum curve to make the comparison easier.

    

    Water is destabilizing the reactants relative to the product and transition state, making it energetically easier for the reactants to reach the transition state and form product. However, the reactants first have to find each other, and then stay together long enough for reaching the transition state. This is an entropic effect which can be estimated by different techniques. We will not do that here, as it is outside the scope of this tutorial.

back to top

Conclusions

We have now found the geometries and energies of the transition state, the reactant state and the prodcut state of the Diels-Alder cycloaddition in water. Table 3 lists the QM/MM energies of these geomtries.

Table 3. Energies of the reactant,
transition state and product geom-
etries in water at the
PM3/GROMOS96 QM/MM level.
The last column lists the energy
differences with respect to the
reactant state.

	E (kJ/mol)	ΔE (kJ/mol)
Reactant	-1202.4	0.0
Trans. St.	-1062.9	139.5
Product	-1271.7	-69.3

back to top

updated 07/09/04