Physikalische Chemie Masterpraktikum Molecular Modelling

Physikalische Chemie
Masterpraktikum
Molecular Modelling
Thermochemistry and reaction mechanisms
- A simple SN2 reaction
Betreuer: Dipl.-Ing. Stephan Kohaut,
Dipl.-Chem. Nathalie Kunkel
Bereiten Sie bitte folgende Themengebiete mit Hilfe des Th2-Manuskriptes vor:
• Hartree-Fock (Kap. 9)
• Basissätze (Kap. 10)
• DFT (Kap. 15 & 18)
• Strukturoptimierung (Kap. 17)
• Schwingungsspektren (Kap. 19)
1
1. Introduction
The exploration of reaction mechanisms and reaction paths that cannot be measured
directly during an experiment has nowadays become a daily routine for chemists to
support their laboratory work. Among all ab-initio techniques Density Functional
Theory (DFT) has become the golden standard due to its wide applicability, ranging
from simple organic molecules to complex coordination compounds involving
transition metals. For the prediction of thermo-chemical values one has to make a
step from a DFT calculation on a single molecule to macroscopic, observable values.
With the help of statistical thermodynamics it is possible to describe an ensemble of
molecules and close the gap between atomic energy units and units that are more
familiar to chemists like for example kJ/mol. The properties of such ensembles are
described by fundamental thermodynamical variables like enthalpy, entropy,
temperature, internal energy etc.. If one finds a way to calculate these values it would
be possible to predict the most important values for synthetic working chemists,
namely the differences in the heat of enthalpy, the energy of activation and the Gibbs
free energy for a chemical reaction. Nevertheless at the beginning one has to
accurately calculate single molecule properties with DFT by using an appropriate
functional. “Which functional should I choose?” is the first question that comes into
ones mind when dealing with DFT calculations [1]. All functionals currently used are
still approximations to the exchange-correlation problem in DFT resulting in the fact
that no functional is accurate for all systems and all wanted properties. This fact
should be covered by the DFT “Jacob’s ladder to heaven” on which every rung
represents a different and more complex approximation, that should have the
features of former rungs with additional wider applicability to new problems [2].
Starting from local density approximation functionals (LDA) that are rather useless
for calculations with chemical accuracy due to strong overbinding tendency to
generalized gradient approximations (GGA) that partially correct this fault, to metaGGA’s and hybrid functionals that are at present the most widely used functionals for
molecular DFT calculations. One should also keep in mind that functionals do not
perform equally well for molecular systems and crystals with their translational
symmetry properties.
1
While the local functional LDA performs well for the prediction of for example lattice
constants in crystals due to their slow varying electron density, it fails for
inhomogeneous systems that are present in molecular structures.
By far the most popular functional is the hybrid B3LYP [3,4] which contains about
20% of exact Hartree-Fock exchange to mimic effects of static correlation:
(1)
Yet today's functionals do not reach chemical accuracies for the prediction of
reaction energies where typical errors are in the range of 3-5 kcal/mol [1].
Why are then some functionals more popular than others? The answer is connected
to the question which system I want to describe or more precisely in which specific
property one is interested! Nevertheless a good starting point are always functionals
that are applicable to lots of different questions, which automatically results in a loss
of accuracy. Although their results can be crude, they often give at least a
qualitatively correct result as a starting point to further calculations. Such functionals
that show wide-range applicability are PBE (GGA) and the already mentioned
B3LYP (hybrid).
2. Thermochemical calculations [2,3]
All calculations assume an ideal gas of non-interacting particles in the rigid rotor,
harmonic normal mode approximation. It is assumed that the first and higher
electronic excitations are not accessible, a fact which could be troublesome in the
case of the existence of low lying electronic excited states. During this exercise the
changes in enthalpy and in the Gibbs free energy together with the energy of
activation will be calculated for a simple SN 2 reaction at room temperature as shown
in equation 2.
(2)
2
Usually chemists are exclusively interested in stationary points (local extrema),
where all forces vanish , whereas the properties of such a point are determined by its
Hessian Eigenvalues. To distinguish between a ground state and a transition state or
saddle point, one has to calculate the mass-weighted Hessian matrix for a given
molecular configuration. Murrel and Laidler defined a transition state as having a
single negative Hessian Eigenvalue corresponding to a negative force constant and
a single negative imaginary normal mode frequency [5]. Transition states are rather
important in terms of chemical reactions because they can be attributed to low
energy paths between two ground states (fig.1) leading to the possibility to calculate
energy of activations and rate constants from transition state theory (fig.2).
Figure 1: An example of a simple PES with two minima (M) connected by a transition
state (TS) and a higher order saddle point (S) [6]
Figure 2: Energy profile as a function of a
reaction coordinate with two minima and a TS
[7]
3
The enthalpy E at a temperature of 298K can be calculated from equation 3 [8]:
(3)
Whereas
can be calculated from the following equation 4 [8]:
(4)
= Difference in energy between educts and reactants at 0K
= Difference in the zero point energy between educts and reactants at 0 K
= Change in the vibrational energy by going from 0K to 298K
= Difference in the rotational energy between products and reactants
= Difference in the translational energy between products and reactants
dN with dN beeing the difference in the number of molecules of reactants
and products
The free Gibbs energy
can be calculated from equation 5 as:
(5)
where
contains all contributions to the total partition function corrected to 298K
(see at ref.9):
(6)
The program will calculate and print all corrections to the thermal energy, the
enthalpy and the Gibbs free energy at the end of every *.out file but one has to take
care and check for different units used (convert them!).
4
Remember that the values corrected to 298K are already containing the zero point
energy (ZPE)! Never add the ZPE except for calculations at 0K, in that case it must
be added manually.
To calculate the reaction enthalpy and the Gibbs energy one simple has to use
equation (7) and (8) by using the corrections to enthalpy and Gibbs free energy and
the printed total energies for every molecule [9]:
(7)
(8)
It is recommended to print a table containing all the calculated total energies,
enthalpies and Gibbs free enthalpies for every species involved in the reaction. How
is the energy of activation calculated then?
3. Instructions
● Optimize the geometries for all educts and products with DFT(B3LYP) and a
6-311++G** basis set [10] and compare them with published results. Use
point groups for every species involved to speed up all calculations.
● What are the geometries for [GS1]-, [GS2]- and [TS1]-? Do they correspond
with published results? To find the geometry [TS]- start with a “good” guess
based on your chemical intuition together with the proper point group.
● Verify that every structure is indeed a ground- or transition state by running a
mass-weighted normal coordinate analysis at the same level of theory.
● Calculate the change in enthalpy, Gibbs free energy and the energy of
activation (for educt -> TS1 -> product) for every reaction step at 298 K and
print a table with all results and a diagram of Gibbs free energy vs. reaction
coordinate
(enthalpies
in
brackets!).
spontaneously?
5
Does
the
reaction
happen
● Calculate the reaction rate with the Eyring–Polanyi equation and compare
with results from the literature (most likely an MP2 calculation).
● If there is time left verify that TS1 connects your reactants and products by
doing an intrinsic reaction coordinate (IRC) calculation
4. Manual
In every student's folder will be an input template file *.inp containing all the basic
keywords for doing calculations with the firefly program [11]:
$CONTRL RUNTYP=optimize maxit=200 exetyp=run nprint=9 $END
$CONTRL SCFTYP=RHF $END
$CONTRL DFTTYP=B3LYP $END
$CONTRL ICHARG=0 MULT=1 d5=.true. icut=11 $END
$SCF DIRSCF=.t. damp=.t. diis=.t. fdiff=.f. soscf=.f. shift=.t. nconv=7 $end
$SYSTEM MWORDS=175 timlim=60000 $END
!$IRC saddle=.t. tsengy=.t. forwrd=.t. npoint=100 stride=0.1 opttol=0.0005 $end
$p2p p2p=.t. dlb=.t. $end
$mp2grd tol1=1d-12 tol2=1d-12 $end
$statpt nstep=300 hssend=.t. OPTTOL=0.000009 $end
$Force nvib=2 $end
To draw a molecule open the chemcraft program and press strg+a (prints atoms) or
strg+f (various molecules ready for usage). Afterwards orientate the molecule in the
corresponding point group by choosing Edit - set point group ( if you are not sure
which point group is the right one you can try the automatic search, but never use
this feature without using your brain!). If the molecule has the correct point group
symmetry go to Tools - Create section for input file - GAMESS US and choose the
basis set 6-311++G**. Press the button copy and insert everything below the
template input. It should be formatted and looking like this:
$Force nvib=2 $end
$CONTRL COORD=UNIQUE $END
$DATA
6
HO
Cnv
8
first atom
8.0
0.0000000000
0.0000000000
0.00000000000
….rest....
By choosing RUNTYP=optimize a geometry optimization will be done which ends if
the RMS of the gradient on all atoms falls below OPTTOOL. Afterwards a frequency
calculation will be performed due to hssend=.t. in $statpt. All results including the
frequency calculations can be visualized with chemcraft by opening the *.out file.
As mentioned before the numbers needed for the calculation of the thermochemical
values are at the end of each *.out file, they can be opened by using a text editor (for
example gedit). To search for TS1 you must first draw a guess, then change the
input keyword from RUNTYP=optimize to RUNTYP=sadpoint and add to $statpt
nstep=300 hess=calc OPTTOL=0.000009 $end. This will perform a normal
coordinate analysis following a geometry search along the largest imaginary normal
mode to locate the TS point. If there is time left, start an IRC calculation to see if your
found TS geometry connects your reactants and products. For doing that change
RUNTYP=optimize to RUNTYP=irc and copy the $HESS group after a frequency
calculation from the file PUNCH to the end of your input file after $END. All results
can be visualized within chemcraft by using the View tab and the proper flags there,
the molecules itself can be saved as a picture file by using File - save image.
5. Sources
[1]
Rappoport ,D.; Crawford ,N.M.R.; Furche , F.; Burke,K. Computational
Inorganic
and Bioinorganic Chemistry ; E. I. Solomon, R. B. King, and R.
A. Scott, Eds.,Wiley, Chichester (2009)
[2]
Perdew,J.P.; Schmidt, K.; AIP Conf. Proc. 577, 1 (2001)
[3]
Becke,A.D. J.Chem.Phys. 98 (1993) 5648-5652
[4]
Lee, C.W. Yang; Parr,R.G. Phys. Rev. B 37 (1988) 785-789
[5]
Murrel,J.N.; Laidler,K.J.; Trans. Faraday Soc. 64 (1968) 371
7
[6]
Wales, D.J.; in Atomic Clusters and Nanoparticles 437-507 (2001)
[7]
http://en.wikipedia.org/wiki/File:Quasi-equilibrium1.jpg
[8]
Foresman,J.B; Frisch,A.; Exploring Chemistry with Electronic Structure
Methods, Sec.Ed. Gaussian Inc. (1996)
[9]
Ochterski,J.W.; Thermochemistry in Gaussian, Gaussian Inc. (2000)
[10]
Hehre, W. J.; Random, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio
Molecular Orbital Theory; Wiley: New York, 1986.
[11]
Granovsky,A.A;Firefly version RC 8.0.0G,
http://classic.chem.msu.su/gran/firefly/index.html
8