Labs 9 and 10_ The use of solvation models and the ONIOM layered by hcj


									Labs 9 and 10: The use of solvation models and the ONIOM layered approach in

In this lab we will consider two techniques that are very useful to model larger systems:
the use of solvation models to mimic systems in solution, and the use of the ONIOM
model to calculate different parts of the molecule at different levels of accuracy.
Unfortunately they are mutually exclusive in the current implementation in Gaussian.
Moreover the solvation model cannot be used in conjunction with semi-empirical
methods. This is a limitation of Gaussian. There is no good reason why this would not be

A. Continuum solvation models in Gaussian03.

The first type of calculation uses an implicit solvation model. It is known as the polarized
continuum model (PCM) and was developed by Prof. Jacopo Tomasi from the University
of Pisa. In this approach the interaction between a solvent and solute is modelled by a so-
called Self-Consistent Reaction Field model. The molecule is placed inside a cavity,
made up from the van der Waals sphere of the atoms in the molecule. Inside the cavity
we have the permitivity of vacuum, outside the sphere we have a dielectricum with a
dielectric constant that mimics the solvent of interest, for example water has   78.39 ,
while ether has   4.335 (in the Gaussian program). In the calculation the dipole of the
solute interacts with the electric field of the solvent, and the electronic structure of the
molecule is determined in the presence of this interaction. In the PCM model the electric
field by the solvent is represented by apparent surface charges on the boundary of the
cavity. This provides a very effective way of solving the equations. The use of a proper
cavity is the most tricky part of the calculation, but it is handled automatically by
Gaussian. Let us consider an example where solvation effects are important. We will
consider various isomers of glycine, both in the gas phase and in water. Consider the
following structures

              H                  O
    O                                      O                      O   -
                                      C                                    O
          C                                                           C
                             H        C
                      H                                                         H
H                                          N                  H       C
                                 H                                        N+
               N                                                  H
    H                                  H
                                               H                      H
                      H                                                     H

        Isomer A                     Isomer B                     Zwitter Ion
Let us first build these three versions of glycine in the gas phase, pre-optimize structures
at the AM1 level, and then run full optimizations at the B3LYP / 6-31G(d) level. For
more accurate structures we should really use a better basis set (e.g. 6-31++G(d,p)), but
this takes too much time for our lab. Our optimization provides gas phase geometries.
Next start from the B3LYP optimized gas phase structures and reoptimize the geometry,
including effects of solvation. Under solvation select the default model (a variety
of PCM) and select water as the solvent. Pull down the menu to see the other
possibilities. The solvents are arranged in order of decreasing polarity (dielectric
constant). Now we can optimize the three versions of glycine, mimicking the aqueous
environment. The energies we get are not really good enough yet, because the basis set is
lacking. For each of the optimized geometries and solvation models (gas phase or water)
run single point energy calculations using the 6-31++G(d,p) basis set. This includes
diffuse functions on all atoms (including H) and polarization functions on hydrogen also.
This is in particularly important for the gas phase calculations. Moreover, on the keyword
line put SCF=Tight. Gaussian prints out a warning, that this is needed with diffuse
functions otherwise. Now we can fill in relative energies for these isomers in the gas
phase and in solution.

Assignments on solvation effects:

1. Finish the above exercise and calculate the relative energetics in KJ/mol for the above
three structures both in vacuum and in water. Now repeat the set of calculations using
ether as a solvent (both structure optimization at the B3LYP / 6-31G(d) level and higher
accuracy single point calculations (B3LYP / 6-31++G(d,p) with SCF=Tight). Comment
on your findings.

2. Also find the transition state for the hydrogen transfer between the structure of Isomer
B and the Zwitter ion. I found it most convenient to first do a QST3 calculation at the
AM1 level, calculating the force constants in every step. Using the optimized AM1
transition state you can then locate the transition state (using QST3) at the B3LYP/6-
31G(d) level, both in the gas phase and in solution (water solvent). Calculate the
frequencies to check you have indeed located the transition state. Once you have found
the transition state calculate the energetics at the higher level (B3LYP // 6-31++G(d,p),
SCF=Tight). Now you can fill in the electronic energies for all of the structures. To make
comparison with experiment we would need to also calculate vibrational frequencies to
include the effects due to zeropoint motion and other thermal effects, but this time around
let’s just say we did …
B. ONIOM calculations in Gaussian03.

Let us first do a warm-up problem to show the use of the ONIOM approach. Consider the
calculation of the carbonyl stretch frequency in propanal, CH3CH2CHO.

                                         Low Level
                                                                            Link atom
          H                  H
                     C                                                H          H
               C                 H
         H                   C                 High level
                             O                                             O

In the (double layered) ONIOM approach we divide the molecule into two pieces. The
region of most interest (containing the carbonyl group here), and the rest of the molecule.
If we are interested in the CO bond length and frequency we can approximate the energy
        E(ONIOM) = E(propanal / low level) +
                E(formaldehyde / high level) – E(formaldehyde / low level).

We can use this compound formula for the energy to optimize the structure and to
calculate the force constants. The positions of all the atoms are optimized for the
complete structure, while the position of the link atom is somewhat arbitrary and uses a
protocol developed in the group of Keji Morokuma, the inventor of ONIOM, at Emory
University. Such a calculation is quite easy to set up in Gaussview. First draw propanal.
Then select Multilayer Oniom model in the calculation set up under the method
card. We now have to assign the atoms to the low and high level regions. This is done
using the layer selection tool under the edit menu. You might first select all
atoms in the low region (use set layer), and apply this selection. Next go on to the
high layer, first initialize by select none and then select atoms by clicking on them
and holding down the shift key. When you are done selecting the atoms use apply to
make the changes permanent. In Gausview the atoms in different layers are rendered
differently, so you can check your selection. Under
view/display/format/molecule you can change the visualization. I prefer tubes
in my low layer, in case of a double layer calculation. At this point we have selected the
layers. In setting up the calculation we now have to specify the type of calculation in each
layer. For our current application we may select the MP2/6-31G(d,p) for the high
accuracy region, and HF/3-21G or AM1 for the low accuracy calculation.
Gaussian will specify the link atom automatically for this simple calculation. It is useful
to look at your input file ( before running the calculation. You will see that
each atom is assigned a level H or L. The carbon atom which is the link between the low
and high accuracy regions has additional “H 1” on the input line. This specifies that in the
high accuracy calculation the atom is replaced by H and it is bonded to atom 1 (the
carbonyl carbon). Your numbering of the atoms may be different, it depends on your
construction of the molecule. To see the atom numbering, use view/labels. This will
show how the atoms are numbered in the input file. If the input file does not create the
proper link atom, I think it is easiest to fix this in the input file directly. There are other
possibilities to do this, which I think are more complicated. You can run the calculation
and visualize frequencies as before. You will see that the carbonyl stretch does include
some motion in your low accuracy region. We can increase the size of the calculation by
selecting also the CH2 group to be part of the high accuracy region. The calculation will
be more expensive, but you expect it to converge to the result we would get if we would
do a complete MP2/6-31G(d,p) calculation. Let us finally extend the model one more
step, and use a three-layer ONIOM model. We use formaldehyde in the high layer, using
MP2 / 6-31G(d,p), then select the CH2 group for the second, medium layer, and will do a
HF / 6-31G(d,p) calculation there, while finally we will use HF/3-21G calculations for
the final CH3 group. This example is a little small for such a complicated approach of
course, but it is easy to envisage how such an approach can be useful. It is definitely
advantageous to make gradual changes in the computational level.

Assignments on ONIOM:

1. The propanal molecule is small enough that we can also perform a full MP2 / 6-
31G(d,p) calculation. List the frequencies corresponding to the frequencies in the high
accuracy region, and list the relevant geometrical parameters. Compare these quantities
for all four (ONIOM) calculations that we did above, and identify the most appropriate
ONIOM calculation.

2. Consider prototype electrophylic aromatic substitution reactions involving attack on
benzene by either the Chlorine cation or the NO2 cation. These reactions involve
intermediates (i.e. minima on the potential energy surface) of the following form

where X = Cl, or X = NO2. We have seen these before. Don’t worry; no more transition
states ... The intermediates have charge 1, and they are spin singlet states.

a) Optimize the structures at the AM1 level and investigate the (Mulliken) charges.
b) For both species substitute either the para or one of the meta hydrogens by a t-butyl
group (one of the predefined groups in Gaussview). The t-butyl group is electron
donating and you can expect that it will stabilize the intermediate if is attached to the
carbon that has a positive charge in the parent cation under a). Let us investigate this by
performing an ONIOM calculation. Select the methyl groups of the t-butyl group as your
low level and the remainder (including the central carbon of the t-butyl group) as the
high-level part of the ONIOM. We can use DFT/6-31G(d) for the high level and AM1 for
the low level. Optimize your 4 structures (meta/para for either Cl or NO2 species),
compare the energies of the para and meta substituted species and comment on the fact if
this agrees with your expectations, based on the charge distributions seen under a). Also
examine and compare the geometrical parameters involving the H-C-X moiety.

c) For graduate students only: Carry out improved single point energy calculations at
the ONIOM optimized geometries. We might use ONIOM with B3LYP/CC-PVTZ for
the high accuracy region and 3-21G for the low accuracy region. Finally replace the t-
butyl group by CH3 and carry out optimizations at the B3LYP/6-31G(d) level, followed
by B3LYP/cc-PVTZ single point calculations (i.e. no ONIOM). Compare the geometries,
and energetics between the methylated and t-butylated species. Also compare the charges
on the most relevant atoms. I would not expect this substitution to have a huge impact, in
particularly not on geometries. If the results differ substantially, it is more likely to be an
error introduced by the ONIOM calculation on the t-butylated species. It is not entirely
trivial to mix the proper computational levels and get accurate results from ONIOM. If
you perform ONIOM calculations in actual research you should properly test things out,
by making comparisons between gradually more accurate approaches.

To top