In this lab we will examine the use of the second derivative of
the energy with respect to nuclear displacements in computational
chemistry. Second derivatives have a number of uses including
identification of a stationairy point as a minimum, a transition
state or something else, the evaluation of the entropy contribution
to Gibbs free energy and the calculation of vibrational spectra.
In the previous two labs we have studied calculations where the
first derivative of the energy with respect to geometry were of primary importance.
In a geometry optimization we look for a minimum in the energy. At
a minimum the first derivative of the energy with respect to any
change in the geometry is zero and the energy will go up
with any change in geometry. At a transition state the first
derivative of the energy with respect to any
change in the geometry is again zero and the energy will go up with any
change in geometry except in one direction where the energy will
go down (the MEP).
The first derivative is the slope of the PES
at a certain point. The second derivative is the curvature
of the PES at a certain point. As we shall see, the second derivative
contains useful information. It also takes longer to calculate than
the first derivative.
In infrared spectroscopy the observed signals are usually
associated with molecular vibrations.
From your physical chemistry classes you may recall that the
energy required (in cm-1) to cause a vibrational mode of a diatomic molecule to
increase in quantum number by one is
E(n→n+1) =
h
2πc
√
k
m
(6.1)
where h is Planck's constant c is the speed of light,
k is the force constant associated with the vibration
and m is the reduced mass. h and c are physical constants,
2 and π are numbers and m depends on the masses of the
two atoms that make up the diatomic. The final variable is the force constant k.
k is simply the second derivative of the energy at the minimum in
the potential energy curve. Thus, if we
know the second derivative of the energy of a diatomic molecule at this minimum
and the masses of the atoms in the molecule then we can
evaluate the energy of the vibration of the molecule.
The process is somewhat more complicated for a polyatomic
molecule because the form of the vibrations (the normal
modes) must also be determined but this can be done.
You may also recall that the intensity of a band
in an infrared spectrum depends on the change in dipole
moment of the molecule as the vibration takes place or, in
an equation
Ii ∝
∂μ
∂Qi
(6.2)
where Qi is the normal mode of the ith vibration.
Once the second derivatives of the energy have been calculated
the ¶m/ ¶Qi can also be calculated
quite easily.
With the energies and intensities of the vibrational
excitations available it is possible to predict
what the IR spectrum of the molecule will look like.
Further Information
For more theoretical background on
infrared spectroscopy look at
http://orgchem.colorado.edu/hndbksupport/irtutor/tutorial.html
or http://physics.schooltool.nl/irspectroscopy/theory.php.
A detailed theoretical description is given in the book
Molecular Vibrations by E.B. Wilson Jr, J.C. Decius and
P.C. Cross, Dover, New York, 1980. A nice description of theory
required to calculate infrared spectra using quantum mechanics
is given in the article by J. Neugebauer, M. Reiher, C. Kind and
B.A. Hess in J. Comput. Chem. vol. 23 pg 895-910 2002.
We choose to calculate the second derivatives of the energy
by selecting the Frequencies option from the
Preset menu on the Main Options window (Figure 6.1). As usual,
some option of the frequencies calculations can be modified
through the Task:Frequencies option of the
Details menu but usually the default parameters will work
well.
Figure 6.2: A SPECTRA window after a frequencies calculation.
Once the second derivatives have been calculated the predicted
IR spectrum can be analyzed using the SPECTRA program from the
SCM menu. If a Frequencies calculation has been run we can
run SPECTRA and a simulated spectrum should appear (figure 6.2).
If the mouse is moved over a peak the energy and intensity of that
peak will appear. It is also possible to visualize what a given vibration
looks like using the SPECTRA program. If the left mouse button is pressed
while the arrow is over a particular peak then a window
(or windows for degenerate vibrations) will appear showing a
movie of the vibration.
Some peaks may be difficult to find because their
intensity is low. Other peaks will have zero intensity by symmetry.
The energies and intensities of these peaks can
be accessed through menus. Each symmetry irreducible representation
that has a mode has its own menu. Clicking on this menu will list all
vibrations of that symmetry and their energies and intensities. Choosing a
vibration from the list will create a red dot on your simulated spectrum.
Chicking on this red dot will produce a movie of the vibration.
Of the other options of SPECTRA, those along the
bottom of the window control the shape and width of the peaks,
the options in the axes menu control the units and overall
form of the spectrum and those in the file menu allow the spectrum
to be saved (as a picture or as pointwise data).
Exercise 6.1
Use a DZP basis set and
the LDA functional for this exercise.
Optimize the geometry of ammonia, NH3. Evaluate the
second derivatives of the energy at the optimized geometry. Compare the length of time
that it took to optimize the geometry of the molecule with the
length of time that it took to calculate the second derivatives. Simulate
the infrared spectrum of NH3 with SPECTRA. What is the general
formula for the number of vibrations of a nonlinear polyatomic molecule with
N atoms? Do you get this number of peaks in your simulated spectrum? Why or why not?
Save a picture of your spectrum in the transmission form with units
cm-1 and peak width 100. Save a second picture of your spectrum in
the absorbtion form with units eV and peak width 0.003. Which peak corresponds
to a scissoring motion of the H atoms? Which peak corresponds to the
symmetric stretching of the N-H bonds?
6.4 Characterization of Minima and Transition States
As was mentioned in section 6.2, a minimum and a transition state
are similar in that the first derivatives with respect to all
changes in geometry are zero for both. Other
situations exist where the first derivatives are zero but we don't
have a minimum or a transition state. It is possible
to ask for a transition state and find a minimum because, in terms of first
derivatives at least, they both look the same. If a geometry
optimization or a transition state search finished
successfully it is possible to say that
the first derivatives are zero but not for sure whether what has been
found is a minimum or a transition state or something else.
Second derivatives can be used
to confirm that what has been found is a minimum, a transition
state or something else. You may recall from calculus that if a
one-dimensional function is minimized the second derivative of
that function will be positive. If the second derivative is
negative then the function is maximized. (If you don't
recall this, take a look at
http://www.ugrad.math.ubc.ca/coursedoc/math100/notes/apps/second-deriv.html,
http://en.wikipedia.org/wiki/Maxima_and_minima
or http://www.karlscalculus.org/calc5_0.html)
Similar ideas hold for functions of more than one dimension like
a PES. Here we have many dimensions. If the gradient of the energy
is zero in all dimensions then the energy is a maximum or a
minimum in each direction. Each direction has a second derivative and
if that second derivative is positive then the energy is minimized
in that direction while if the second derivative is negative then the
energy is maximized in that direction. Therefore, at a local minimum
the second derivative is positive for all directions and at
a transition state the second derivative is positive in all
directions except one.
Figure 6.3: The SPECTRA window of a transition state. Note
the one negative intensity frequency.
Evaluating the second derivatives is clearly a good test
to see if the stationairy point (place where all first derivatives are zero)
that has been found is a minimum (all second derivatives > 0) a
transition state (exactly one second derivative < 0) or
something else ( > 1 second derivative < 0).
In the SPECTRA program vibrations with second derivative
< 0 are easy to see because they have energy and intensity
< 0. The movie of the negative vibration is particularly
informative because it shows how the geometry of the molecule
changes as it crosses the TS state of the reaction.
Note that in reality the frequency of a vibration that has
a negative second derivative is actually imaginary because
E(n→ n+1) √ k. Imaginary numbers are
hard to plot however so negative works better.
Exercise 6.2
In the last lab we considered the reaction
HCN → CNH. We calculated the reactants, the products
and the transition state. Take your optimized geometries of
the reactants, products and TS and run a frequencies calculations
for each to confirm that the first two are indeed minima
and the third is indeed a transition state. List the vibrational frequencies of
the reactant, product and TS. Briefly describe the geometry change
occurring at the transition state of this reaction.
Many of the equations of this section come from the theory of
statistical mechanics. We will not go over any of the details of how
the equations come about. If you are interested in the details of
this theory see your physical chemistry text books.
In the previous lab we calculated the energies of reactants, products
and transistion states of reactions and used these results to approximate
reaction enthalpies and enthalpy barriers. Enthalpies and free energies
are macroscopic quantities that can be directly measured. We calculate
the properties of a single molecule. Using statistical mechanics the properties
of one molecule can be extrapolated to the macroscopic frame. Now that we have
the energy second derivatives we have all the information required to calculate enthalpies.
What we have been calculating so far is the electronic contribution to
the enthalpy. The enthalpy can be written as
H = Eelec + EZPE + Eint +nRT
(6.3)
where Eelec is the electronic contribution that we already have, EZPE
is the vibration zero-point energy, Eint is the internal energy and
nRT is the energy of the system an constant pressure as is required
for the enthalpy.
EZPE depends on the vibrational frequencies of the molecule. Once we
have calculated the second derivatives we have these frequencies and so can
evaluate EZPE
ADF calculates EZPE whenever a frequencies run is performed. It can be found
in the detailed output. If the OUTPUT program is run the zero-point
energy can be found through the Other Properties:Zero-point Energy menu
item (see figure 6.4).
Eint is made up of thermal energy divided into contributions from
translational, rotational, vibrational and electronic motion. We will always assume
that the electronic thermal energy is zero. The evaluation of the translational thermal
energy requires that we know the mass of the molecule. The rotational
thermal energy can be evaluated if we know the masses of the atoms in the molecule
and the geometry of the molecule. The vibrational thermal energy can be calculated
if we know the masses of the atoms in the molecule and the energies of all its
vibrations. All of this information is available once a frequencies calculation
has been completed. It can be found in the detailed output by selecting
the Other Properties:Statistical Thermal Analysis menu option
of the OUTPUT program. The translational, rotational and vibrational
internal energies are listed in kCal/mol. It is important to note that
the EZPE is included in the vibrational internal energy.
Figure 6.5: Internal energy contributions.
The thermal contributions to the internal energy are obviously temperature
dependent. By default the temperature is chosen to be 298.15 K. You
can choose a different temperature or a range of temperatures through the
Properties:Thermodynamics option of the ADFINPUT program.
Exercise 6.3
Make up a table of the thermodynamic data of the reactants,
products and transition state of the HCN→CNH reaction at 298.15 K.
Your table should include a row for Eelec, EZPE,
Eint, nRT and the enthalpy H. Don't forget that ADF includes EZPE
in the vibrational part of Eint so you should be ensure that you
do not count it twice. From the data in your table,
evaluate ΔH the enthalpy of the reaction and ΔHf the
enthalpy barrier of the reaction. Compare your ΔH values with the
approximate energies you evaluated last lab using just Eelec. You can use any units
for your analysis but give your final ΔH values in kJ/mol.
The enthalpies of a reaction are important but more important still
are the Gibbs free energy and free energy barrier of the reaction, ΔG
and ΔGf. With these quantities in hand it is possible to
evaluate the important chemical quantities the equilibrium constant and
the rate constant since
Keq ∝ e-[(ΔG)/kT]
(6.4)
krate ∝ e[(ΔGf)/kT]
(6.5)
This equation for the rate constant is more correct than the one given
in the previous lab.
We know from thermodynamics that
G = H - TS
(6.6)
We have already gone through how to calculate H. S is similar to the internal
energy in that it has contributions from translation, rotation and vibration.
The information required to calculate these entropy contributions is
the same as that required to obtain the internal energy contributions. Thus,
it is not surprising that the contributions to the entropy of the molecule (in cal/mol K)
are listed in the detailed output in the same place as the contributions
to the internal energy.
(The Other Properties:Statistical Thermal Analysis menu option
of the OUTPUT program( see figure 6.6.)
Exercise 6.4
Make another table where you list the entropies of the
reactants, products and transitions state of the reaction HCN → CNH.
Add to your table the Gibbs free energy of each species. Evaluate ΔG
and ΔGf of the reaction. Compare ΔG and ΔGf
with the values you obtained for ΔE and ΔEf in the last lab.
Could you have predicted ahead of time the qualitative result of this comparison? Why?
Evaluate ΔEelec, ΔH and ΔG for the reaction
H2O + BH3 → H2O-BH3. Note that we are only interested
in the overall reaction energy. We are not interested in the barrier or
transition state. Use the PBE functional and a TZP basis. The calculations
involving BH3-H2O might take a little
while (10 to 20 minutes) so it is a good idea to run them on 2 processors or
4 if you are running outside of normal lab time. How do ΔEelec, ΔH and ΔG
compare with each other? Why the change from the HCN case?
