QM/MM tutorial

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

Introduction

With a good guess of what the transition state should look like, it was rather straightforward to find the transition state geometry. We will now use the more systematic Linear Transit approach to do the same. In the Linear Transit a coordinate is choosen along which the reactants are transformed into product. This so called reaction coordinate is varied linearly while all other degrees of freedom are minimized. Choosing such reaction coordinate requires some intuition and understanding of the process studied, but is in general easier to chose than a reasonable guess geometry. The concept of the reaction coordinate is best explained by an example. In case of a Diels-Alder cyclo-addition, a good reaction coordinate would be the distance between the two atom pairs that are forming the two new bonds upon reaction (figure 3).

Here we will use the gromacs QM/MM features to perform a Linear Transit calculation of the Diels-Alder cyclo-addition reaction (figure 1). The -R group, which was missing in the x-ray model and ignored in the previous part of the tutorial, will now be taken into account. Because this group is unlikely to have a large effect on the reaction, we will describe it at the MM level in our model.

QM/MM subdivision

Figure 4 shows how we split up our system in a QM and MM part. The QM part consists of the same atoms we were using in the previous vacuum calculations and is described at the semi-empirical PM3 level of theory. The remainder of the system, consisting of the tail part (-R, figure 1) is modeled with the GROMOS96 forcefield.

To make use of the QM/MM functionality in Gromacs, we have to

• introduce link atoms at the QM/MM boundary, if needed;
• specify which atoms are to be treated at a QM level;
• specify the QM level, basisset, type of QM/MM interface and so on.

We also need to know how to do a Linear-Transit in Gromacs:

• setting up a Linear-Transit calculation;
• performing a Linear-Transit calculation.

adding link atoms

At the bond that connects the QM and MM subsystems we introduce a link atom. In Gromacs we make use of a special atomtype, called LA. This atomtype is treated as a hydrogen atom in the QM calculation, and as a dummy atom in the forcefield calculation. The link atoms, if any, are part of the system, but have no interaction with any other atom, except that the QM force working on it is distributed over the two atoms of the bond. In the topology the link atom (LA), therefore, is defined as a dummy atom:

[ dummies2 ]
LA QMatom MMatom 1 0.65

specifying the QM atoms

Once we have decided which atoms should be treated by a QM method, we add these atoms, including the link atoms, if any, to the index file. We can either use the make_ndx program, or hack the atoms into the index.ndx file ourselves. The index file we will use in this tutorial is found here: index.ndx. It is possible to constrain the bonds in the QM subsystem along. It is also possible not to constrain them, while the bonds in the MM subsystem are. This is essential for instance if the QM atoms are supposed to undergo bond-breaking/formation reactions. In this case, Gromacs' bondtype 5 is used for the bonds in the QM subsystem:

[ bonds ]
QMatom1 QMatom2 5
QMatom2 QMatom3 5

specifying the QM/MM simulation parameters

The last thing we need to do to setup gromacs for performing QM/MM calculations is to specify what level of QM theory gromacs has to use for the QM subsystem, what QM/MM interface to use, what multiplicity the QM subsystem has, and so on. All these things are defined in the mdp file. The following option lines need to be included for a QM/MM run:

QMMM = yes
QMMM-grps = QMatoms
QMmethod = RHF
QMbasis = 3-21G*
QMMMscheme = ONIOM
QMcharge = 0
QMmult = 1

setting up a Linear-Transit calculation

In this paragraph we describe how to do linear transits with gromacs. Note that Linear Transist calculations have nothing to do with QM/MM and can be done with an MM only description as well. We will refer to the reaction coordinate defined in figure 3.

We want to constrain the distance between the centers of the atompairs involved in the reaction and minimize all other degrees of freedom in the system. The way we can impose our constraint in gromacs, is to put a dummy atom (atomtype XX), with no interaction whatsoever with any other atom of the system, exactly in the middle of the atompairs (figure 3) and constrain the dummy-dummy distance.

[ dummies2 ]
dummy1 atom1 atom2 1 0.5
dummy2 atom3 atom4 1 0.5

[ constraints ]
dummy1 dummy2 2 1.6

The dummy is constructed every step of the simulation/optimization, so that it is always exaclty in the middle between the atompair. Note, in the current version of gromacs, constraints between dummies are not allowed yet, so we will use a little trick here. The trick is explained here. For now, it is not important to know the details of this trick. Furthermore, already in the next release of gromacs it will be possible to apply the dummy-dummy constraints we need here.

performing a Linear-Transit calculation

Now that we know how to constrain the reaction coordinate, we are going to perform minimizations at different constraint lengths. We create different subdirectories, one for every Linear-Transit point. We perform the minimizations in these subdirectories and use the result as input for the minimization in the next subdirectory. After all minimizations along the reaction coordinate have been done, we collect the energies from the subdirectories and plot them as a function of the reaction coordinate.

The procedure is straightforward, but tedious. Therefore, we will make use of scripts. With the first script we create the subdirectories and place the topology files with different constraint legths in them. The second script then goes into the subdirs one by one, runs grompp en mdrun. Finally a third script runs g_energy, which retrieves the energies from the outputfile ener.edr and collects the energies as a function of the constraint lenght. The scripts are available for download:

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

Using these scripts we are going to perform a Linear Transit calculation on the Diels-Alder reaction in figure 1 in vacuo, along the reaction coordinate shown in figure 3.

back to top

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

• The starting structure for this calculation is the transition-state analogue form the x-ray model. We modify this structure first, such that it becomes a good model structure for the reaction. In the previous step we learned how to use molden. We will use it again to change the atomnames, saturate the complex with hydrogen atoms, and add the aliphatic tail. We can skip the first three steps here, as a prepared coordinate, topology and index file are available already. Preparing these files is not difficult at all, but can consume quite some time the first time.
  • import the pdb file of the analogue we have created in step one of the previous part

    localhost:~>molden -P analogue.pdb

  • modify the atomnames in the Z-matrix editor.

  • use the "add line" button to add the hydrogens and the tail part. Note, the tail atoms, except the atom, next to the nitrogen atom should not have Hydrogen atoms, as this part is going to be described with the Gromos96 which employs so-called united-atoms, i.e., the hydrogens are summed into the carbon atoms. The result should look like:

    

  • save the molecule as confin.crd in cartesian xyz format, using the "Cartesian" button.

  • open the confin.crd in an text editor, such as emacs and re-write the confin.crd such that each atom has a unique name. See for example confin.gro.

• We also need to write a QM/MM topology file for this compound. See for example topol_A.itp. Have a look at this file befroe we proceed. For the linear transit we will need a new topology with a different constraint distance for each point of the Linear transit. In our case the scripts mentioned above will take care of this, which is why we can skip this part.

• Furthermore we need an index file, with entries for the atoms that are supposed to be part of the MM and QM subsystems. We can use make_ndx for that, but we will skip this part as well, and use the prepared index file instead. index

• Finally, we create an mdp file with the options for the QM/MM linear transit calculation. We need to specify that we want to energy minimize, how we want to minimize, that we want to do QM/MM, what QM method we want to use. The mdp file we are going to use is LT.mdp.

• With the input files and the scripts at hand, we can now perfom the actual linear transit calculation. We will let the reaction coordinate constraint vary between 0.12 and 0.5 nm in 200 steps.
  • create a new subdirectory for the Linear Transit in vacuo

    localhost:~>mkdir LTvacuo

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

  • 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.12 to 0.5 nm.

    localhost:~>./create_tops.scr 0.12 0.5 200

  • execute the run.scr scripts to run the minimizations in the different subdirectories. 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.12 0.5 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.12 0.5 200 > eqmmm.xvg

• Let's have a look at the results
  • Import eqmmm.xvg in xmgr to see the energy as a function of the reaction coordinate.

    localhost:~>xmgr eqmmm.xvg

    

    The maximum is at 0.2093 nm (-149.924 kJ/mol). This corresponds to point 47. Let's have a look at that structure. We simply go into the directory 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

    

    It looks verly much like the one we optimized in the previous step. You can check if it is a TS by perfoming a frequencies caclulaiton on this structure, using gaussian, just like we did before.

    The minima are at 0.4468 nm (-342.826kJ/mol) and 0.1485 (-388.276 kJ/mol) and correspond to the reactant and product state respectively, and are shown here, along with the transition state again.

    

    Figure 5. Reactant (a), Transition State (b), and Product (c) geometries in vacuo, found with the Linear Transit method. The
    energies of these structures are listed in table 2.

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 vacuo. Table 2 lists the potential energies of these geomtries. A direct comparison with the energies and energy differences found using the optimization routines of gaussian98 (tabel 1) in not valid, because in the current computations the aliphatic tail (-R, figure 1) was taken into account explicitly, while it was ignored in the previous computations.

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

	E (kJ/mol)	ΔE (kJ/mol)
Reactant	-342.826	0.0
Trans. St.	-149.924	192.902
Product	-388.276	-45.45

back to top

updated 28/07/04