# QM/MM tutorial

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

# III. Optimization of product, reactant and transition state geometries in water, using Linear Transit in gromacs

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.

Figure 6. Division of the system in a QM subsystem and an MM subsystem. The QM subsystem is described at the semi-empirical QM level,
while the remainder of the system, consisting of the reactants' aliphatic tail, the water molecules and the Na+ ion, is modeled with the GROMOS96
forcefield.

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

The starting structure for this calculation is the transition-state analogue form the x-ray model, we have modified in the previous part of this tutorial. Then we place this structure in the center of a periodic box and fill that box with 2601 SPC water molecules and 1 Na+ ion. The total system is equilibrated, before the Linear transit computation is performed. As before we will skip this tedious procedure, which has nothing to do with QM/MM, and use the results instead (confin.gro). The steps we took in the equilibration process are described here.

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

The create_tops.scr and the run.scr script are identical to the ones we used before, but the get_energies.scr script is slightly different. In vacuum we could simple take the total potential energy. In the water, and later in the protein, we don't want to know the potential energy of the complete system, but rather the internal energy of the reactans plus the contributions from the interaction of the reactants with the surroundings. So what we want is E(reactants)+E(reactant-solvent)+E(reactants-NA+).

• 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.

# Conclusions

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

updated 07/09/04