1. science
2. mathematics
3. geometry

Chemistry 460
Spring 2015
Dr. Jean M. Standard
March 30, 2015
QUANTUM CHEMISTRY PROJECT 3: ATOMIC AND MOLECULAR STRUCTURE
OUTLINE
In this project, you will carry out quantum mechanical calculations of atomic and molecular structure. In addition to
structure and properties of atoms and molecules, some of the things you will explore are the effects of level of
theory and basis set size on the results. The computer programs that you will use for these simulations include a
state-of-the-art software package called Gaussian for the quantum mechanical calculations and a software package
called Avogadro for visualization.
You will be able to access Gaussian and Avogadro in the computer lab located in Julian Hall Room 216. Please
note that the computers in the JH 216 are usually available at any time during the day; however, they are
occasionally reserved for classes as noted on the schedule posted on the door. In the event that the room is busy
during a time you want to use the Macintoshes, you may use the Macintosh in Julian 222 (outside my office).
You may use any of the Macintosh computers in JH 216 to access the Avogadro/Gaussian programs (instructions
are provided later in this handout). In fact, you can use one Macintosh for one part of the project, and another
Macintosh for another part of the project since none of the data is stored on the Macs. All the calculations actually
will be carried out on Linux workstations in SLB 220 that you will connect to from the Macintoshes via the internet.
DUE DATE
The project is due no later than WEDNESDAY, APRIL 15.
GRADING
Project 3 is worth 50 points and consists of three parts. Part A is worth 15 points, Part B is worth 20 points, and Part
C is worth 15 points.
2
OVERVIEW: QUANTUM CHEMISTRY METHODS
Quantum Mechanics of Many-electron Systems
Because the Schrödinger equation cannot be solved exactly for atoms and molecules containing more than one
electron, the goal in most cases is to obtain the most accurate approximate solution possible within the constraints of
time and available computational resources.
Improving the Accuracy of the Results
There are two primary ways in which the accuracy of the approximate solution may be improved. One approach is
to increase the size of the basis set so that it is closer to the limit of a complete set. Since any arbitrary function
(including the true wavefunction) may be represented exactly using a complete set of functions, increasing the
number of basis functions generally leads to an improved representation of the wavefunction. Another approach to
improve the accuracy of the approximate solution is to increase the level of approximation in the methodology (or
level of theory, as it is commonly stated) employed to obtain the approximate energy and wavefuntion. The two
approaches to improving the accuracy of the approximate solution often are utilized in combination, as illustrated in
Figure 1.
Full CI
y
ur
ac
QCISD
cc
MP4
A
Level of Theory
CCSD(T)
MP2
DFT
HF
Minimal
DZ
TZ
QZ
Basis Set Size
Figure 1. The two key factors that control the accuracy of quantum chemical calculations,
level of theory and basis set size.
Levels of Theory
The lowest level of ab initio quantum theory for atomic and molecular systems is the Hartree-Fock model. This
method is often denoted HF, more specifically RHF for closed shell systems (Restricted Hartree-Fock) and UHF for
open shell systems (Unrestricted Hartree-Fock). This method is characterized by approxmations that lead to the
neglect of electron correlation; that is, electron-electron repulsion is treated only in an average sense.
To improve upon the Hartree-Fock model, one approach for the inclusion of electron correlation is the use of
perturbation theory, which determines a correction to the Hartree-Fock energy. The most common form of
perturbative correction is known as MØller-Plesset Second-Order Perturbation Theory, or MP2. Higher orders of
perturbation theory, such as MP4, also may be used. The MP2 method is the least expensive ab initio method for
the correction of the Hartree-Fock model and inclusion of a significant portion of the electron correlation in the
system. However, the use of the MP2 or MP4 methods comes with additional cost in computing time.
3
Beyond perturbative methods such as MP2 and MP4, the primary methods to improve the approximate solution of
the Schrödinger equation include Coupled-Cluster (CC) theory and Configuration Interaction (CI) methods. Both
of these types of methods move beyond adding corrections to the Hartree-Fock solution and include electron
correlation effects directly into the wavefunction. The wavefunction of the Hartree-Fock method is based upon a
single Slater Determinant. The wavefunctions of the Coupled Cluster and Configuration Interaction methods are
based upon multiple Slater Determinants. There are many varieties of Coupled Cluster theory and Configuration
Interaction; they differ by how many additional Slater Determinants are included in the calculation, which are based
upon excitations from a single reference Slater Determinant (or configuration). The most common Coupled Cluster
methods include CCSD, which adds single and double excitations from the reference, and CCSD(T), which adds
single, double, and a selection of triple excitations. Similary, for Configuration Interaction, common methods are
CISD and CISD(T), which refers to methods which include single, double, and in the latter case, a selection of triple
excitations. Also used frequently are quadratically convergent related methods, QCISD and QCISD(T). In all these
cases, the accuracy improves as the level of excitation increases; hence CISD(T) would be expected to be more
accurate than CISD. The computational cost for the CC and CI methods is prohibitive; thus, these methods can only
be applied to small molecules (generally, systems with fewer than 10 atoms).
Basis Sets
The basis set refers to the collection of functions that is used to represent the approximate atomic or molecular
wavefunction. For molecular systems, the LCAO-MO (linear combination of atomic orbitals to form molecular
orbitals) approximation is generally employed; thus, the basis set consist of a set of atomic functions used to
construct the molecular orbitals. For atoms, the atomic functions are used directly.
The smallest set of atomic functions that can be used in an atomic or molecular calculation is referred to as a
minimal basis set. This set is comprised of atomic functions for all the atomic orbitals occupied in the separated
atoms, and must include all orbitals of a given angular momentum even if the subshell is only partially occupied.
For example, the minimal basis set for carbon atom with electron configuration 1s22s22p2 includes 1s and 2s atomic
functions, along with all three 2p-type atomic functions (2px, 2py, and 2pz), even though the 2p subshell is only
partially filled.
In order to improve the approximate solution of the Schrödinger equation, the basis set size may be increased in
three ways: (1) including more functions with the same angular momentum as those in the minimal basis set but
with different spatial extent; (2) adding polarization functions to better account for the shifts in electron
distributions when atoms are involved in chemical bonding; and (3) adding diffuse functions to account for
electron density far from the atomic nuclei which may be involved in weak interactions such as hydrogen bonding or
van der Waals interactions. Often, the first two approaches for increasing the basis set are employed together, and in
some cases, all three approaches are employed.
If the number of functions employed in the minimal basis set is doubled, we refer to this as a double zeta (DZ) basis
set. If the number of functions in the minimal set is increased by 3×, 4×, or 5×, the basis sets are referred to as triple
zeta (TZ), quadruple zeta (QZ), and quintuple zeta (5Z), respectively. A common double-zeta quality basis set is
6-31G and a common triple-zeta quality basis set is 6-311G; these basis sets are referred to as Pople-style basis sets.
Other common basis sets are referred to as correlation-consistent basis sets; examples include cc-pVDZ, cc-pVTZ,
cc-pVQZ, and cc-pV5X.
Polarization functions in Pople-style basis sets are denoted by including the angular momentum of the basis
functions being added at the end of the basis set name. Examples include 6-31G(d) and 6-311G(d,p). The
correlation-consistent basis functions automatically include polarization functions, and the number and highest
angular momentum of the polarization functions increase as the basis set increases from DZ to 5Z and higher.
Diffuse functions in Pople-style basis sets are denoted by adding '+' after the G in the basis set name. Examples
include 6-31+G(d) and 6-311++G(d,p). The inclusion of diffuse functions in correlation-consistent basis sets is
denoted by adding 'aug' to the name. Examples include aug-cc-pVDZ and aug-cc-pV5Z.
4
USING THE COMPUTERS AND THE AVOGADRO SOFTWARE PACKAGE
This section contains instructions for logging in to the Macintosh computers in JH 216 and accessing the Linux
computers on which the calculations will be performed. In addition, it provides instructions on how to start the
Avogadro software package. The use of the Gaussian software package will be described later.
Logging In
To log in to one of the Macintosh computers in JH 216, enter your ISU ULID and password.
Starting X11
In order to access the Linux workstations to run the Avogadro/Gaussian software packages, you must run a software
package on the Mac called X11. To start the X11 application, go to the Applications folder on the Mac hard drive
and click on the "X" icon. A single white window should appear on the screen; this is called an "xterm" window, or
simply a terminal window.
Logging on to the Linux Computers
You may use any of the six available Linux computers for the project. To log on to a Linux computer, type the
following command in the xterm window:
ssh
–Y
[email protected]
Here, "ssh" stands for secure shell; this application provides a secure connection to the remote Linux workstation.
Also, hostname is the name of a Linux computer [select one of: frodo, samwise, gandalf, aragorn, gimli, or legolas].
The password is "erwin". Enter the password as prompted. Once you have done this, you are connected to the
Linux computer.
Starting the Avogadro Program
To start the Avogadro program, type "avogadro" in the xterm window and hit the return key. A window for the
Avogadro software package should appear on your screen. At this point, you are ready to begin the project.
5
PART A: NOBLE GAS IONIZATION POTENTIALS (15 points)
Introduction
In this part of the project, you will use quantum mechanical calculations to predict the first ionization potentials for
the series of noble gas atoms, investigate trends, and compare the results to experiment. For an atom A, the
ionization potential may be defined as the energy required to remove an electron in the process
A
→
A + + e− .
(1)
In a many-electron atom, Koopmans' theorem states that the atomic orbital energy of the highest occupied orbital
ε HO can be related approximately to
€ the atom's first ionization potential,
IPK = −ε HO ,
€
(2)
where IPK is the ionization potential (from Koopmans' theorem). Koopmans' theorem gives an approximation to
the ionization potential because it assumes that the orbital from which the electron is removed has the same shape
and energy in the neutral species and the cation; however, some relaxation of the orbital usually occurs in the real
system.
Alternately, the ionization potential may be determined from calculations of the total electronic energies of the
neutral atom and cation, Ecation and Eneutral , respectively as
IP = Ecation − Eneutral ,
(3)
where IP is the ionization potential. In this case, the relaxation of the cation orbital is included in the calculation.
Procedure
1.) The first atom that you will perform a calculation on using the Avogadro/Gaussian software packages is helium.
Following the instructions above to start the Avogadro software package should lead to a window similar to the
one shown in Figure 1 appearing on your screen.
Figure 1. The main Avogadro window.
6
The Avogadro window that you see may have a black background. If you would like to change the background
color, you may do so by selecting "View → Set Background Color." In addition, if the "Settings" area on the
left side of the window does not appear, you may toggle it on by clicking on the "Tool Settings" button.
To construct the He atom, click on the Drawing Tool,
. Set the Element to Other and then select He from
the periodic table, Bond Order to Single, and make sure that Adjust Hydrogens is clicked off. Then, click with
the left mouse button in the Avogadro window and a He atom should be displayed in a manner similar to that
shown in Figure 2.
Figure 2. A helium atom in the drawing window.
The next step is to perform a quantum mechanical calculation of the He atom using the Gaussian software
package. The total electronic energy and atomic orbital energies of the atom will be determined. From the
Avogadro menus at the top of the window, select "Extensions → Gaussian". You should see a popup window
similar to that shown in Figure 3.
Make sure that the following settings are entered for the calculation:
Calculation: Single Point Energy
Theory: RHF
Charge: 0
Processors: 4
Basis: 6-31G(d)
Multiplicity: 1
7
Figure 3. Setting up a Gaussian calculation for a noble gas atom.
Once the settings are correct, click "Compute". You will be prompted to enter a filename. Use a filename
something like 'he-XY.com', where 'XY' corresponds to your initials. After the filename is entered, click
"Save". The calculation will begin and you should see a popup window appear that says "Running Gaussian
calculation...". When the popup window disappears, your calculation of the He atom should be complete.
[Note: this calculation is very fast; sometimes you don't see the popup window. If you don't see the popup
window after a minute or two, go on to the next paragraph to check the results file.]
To retrieve the computed value of the highest occupied orbital energy, close Avogadro by selecting "File →
Quit". You may be asked to save changes, but you can safely discard any changes because the file you need is
already saved. Next, in the terminal window, type "gedit he-XY.log" (where 'XY' corresponds to your
initials). You will see a window open on the screen similar to that shown in Figure 4. The file you are viewing
contains the results from the Gaussian calculation that you ran for helium.
8
Figure 4. Results file from a Gaussian calculation of the He atom.
To find the computed value of the orbital energies, search the file for the string 'Alpha occ. eigenvalues' by
selecting "Search → Find" from the menu. The listing you find gives the orbital energies of the occupied and
virtual (unoccupied) atomic orbitals in atomic units (hartrees). For example, for He, the listing given is:
Alpha
occ. eigenvalues
--
–0.91413
Alpha virt. eigenvalues
--
1.39986
For He, there is only one occupied orbital (1s) and its computed energy is –0.91413 hartrees. For each atom,
you should record the energy of the highest occupied orbital.
When you are finished, select "File → Quit" to close the results file.
2.) Carry out similar quantum mechanical calculations for the Ne, Ar, and Kr atoms using Avogadro/Gaussian and
the same procedure described in step #1. You should also perform a similar calculation for Xe; however, in this
case, the basis set selected should be LANL2DZ (all the other settings should be the same).
Once the calculation of each atom completes, find the highest occupied orbital energy in the results file and
record its value.
3.) The next step is to carry out some additional calculations on neon. The first thing you need is the total
electronic energy of Ne from the HF/6-31G(d) calculation you performed in step #2. Use the 'gedit' command
to open the log file for this calculation again.
9
To obtain the total electronic energy, scroll to the bottom of the file, and look for a block of text similar to the
following:
1\1\GINC-FRODO\SP\RHF\6-31G(d)\Ne1\STANDARD\22-Mar-2015\0\\#n RHF/6-31G(d)
SP\\Title\\0,1\Ne,0,-3.3764,1.89097,0.\\Version=AS64L-G09RevD.01\State=1A1G\HF=-128.4744065\RMSD=2.635e-12\Dipole=0.,0.,0.\Quadrupole=0.,0.,0.,
0.,0.,0.\PG=OH [O(Ne1)]\\@
The total electronic energy is listed in bold above (HF=xxxxx) and is given in hartrees. Record this value.
4.) Now repeat the calculation of the neon atom using different levels of theory and basis sets. You will probably
want to do this by building new copies of the neon atom with Avogadro and changing the settings for the
Gaussian calculation as appropriate. It is also helpful to save each calculation with a different filename in case
you need to refer back to them. Once you have selected "Extensions → Gaussian", use the following settings
for the calculations:
Calculation: Single Point Energy
Theory: RHF
Charge: 0
Processors: 4
Basis: 6-31G(d)
Multiplicity: 1
You should see on the second line in the white input box the following text:
#n RHF/6-31G(d) SP
Edit the line so that it has the following format:
#n RHF/aug-cc-pVDZ SP
Note that here you are changing the basis set from 6-31G(d) to aug-cc-pVDZ. In addition, you should add a
line at the beginning of the white input box that has the form:
%mem=1GB
Once you have made all the changes in the settings, click "Compute" to submit the calculation. When the
calculation is complete, follow the instructions given in step #3 to retrieve and record the total electronic energy
from the results file. Make sure to note the level of theory and basis set of the calculation.
Next, repeat the calculation for additional times, changing the level of theory and basis set to those shown in the
table below. Open the results file for each calculation using the 'gedit' command and record the total electronic
energy of neon for each calculation. When you are finished, you should have results for the following six
calculations of the neon atom:
Run
Theory
Basis Set
A (original, step #1)
RHF
6-31G(d)
B (above)
RHF
aug-cc-pVDZ
C
MP2
aug-cc-pVDZ
D
CCSD(T)
aug-cc-pVDZ
E
CCSD(T)
aug-cc-pVTZ
F
CCSD(T)
aug-cc-pVQZ
10
5.) The next step is to repeat the six calculations for the neon cation. You should build each neon cation using
Avogadro as you did in step #1, but build them as if you were making a neutral atom because the charge of +1
will be set before the calculations are performed. From the Avogadro menus at the top of the window, select
"Extensions → Gaussian". You should see a popup window similar to that shown in Figure 5.
Figure 5. Setting up a Gaussian calculation for neon cation.
For the neon cation calculations, make sure that the following settings are selected from the pull-down menus:
Calculation: Single Point Energy
Theory: RHF
Charge: 1
Processors: 4
Basis: 6-31G(d)
Multiplicity: 2
You should see on the second line in the white input box the following text:
#n RHF/6-31G(d) SP
Edit the line so that it has the following format:
#n UHF/6-31G(d) SP
Note that here you are changing the level of theory from RHF (Hartree-Fock theory for closed shells) to UHF
(Hartree-Fock theory for open shells). In addition, you should add a line at the beginning of the white input box
that has the form:
%mem=1GB
Once you have made all the changes in the settings, click "Compute" to submit the calculation. When the
calculation is complete, follow the instructions given in step #3 to retrieve and record the total electronic energy
from the results file. Make sure to note the corresponding level of theory and basis set.
11
Next, repeat the calculation of the neon cation, changing the level of theory and basis set, in order to complete
five additional calculations as shown below. Open the results file with 'gedit' and record the total electronic
energy of the neon cation for each calculation. When you are finished, you should have results for the
following calculations of the neon cation:
Run
Theory
Basis Set
A (above)
UHF
6-31G(d)
B
UHF
aug-cc-pVDZ
C
MP2
aug-cc-pVDZ
D
CCSD(T)
aug-cc-pVDZ
E
CCSD(T)
aug-cc-pVTZ
F
CCSD(T)
aug-cc-pVQZ
Results, Analysis, and Discussion
1.) Use Koopmans' theorem to predict the first ionization potential for each noble gas atom using Equation (2) and
the RHF/6-31G(d) calculations that you completed in step #1. Report your results in both hartrees and eV.
2.) Obtain literature values of the experimental first ionization potentials of He, Ne, Ar, Kr, and Xe. Remember to
cite your source; you should also record the ionization potential for radon because you will need it in #5.
Compare the experimental values to the calculated values. What are the percent errors?
3.) Construct a graph in which your calculated ionization potentials from Koopmans' theorem are plotted on the yaxis and the atomic number of each atom is plotted on the x-axis. Include this plot in your project.
Discuss qualitatively any trends observed in first ionization potentials as a function of atomic number. Provide
physical reasons for the observed behavior.
4.) Construct another graph in which your calculated ionization potentials from Koopmans' theorem are plotted on
the y-axis and a literature value of the atomic radius of each atom is plotted on the x-axis. Again, remember to
cite your source for the atomic radii; you should also record the atomic radius for radon as you will need it in
#5. Include this plot in your project.
Discuss qualitatively any trends observed in first ionization potentials as a function of atomic radius. Is the
trend linear? Fit the observed trend to a polynomial or power series and give the equation of the trendline.
5.) Obtain the atomic radius of radon from the literature. Use the best-fit equation obtained in #4 along with the
atomic radius of radon to obtain a prediction for the ionization potential of radon. Compare to a literature value
and discuss the agreement.
6.) Use the results from steps #4 and #5 of the procedure to determine the ionization potential of neon at different
levels of theory. Use Equation (3) to perform the calculations. Tabulate your results along with percent errors
compared to the experimental value.
7.) Discuss your findings with respect to variations in level of theory with a fixed basis set of aug-cc-pVDZ. Are
there trends in the percent errors as the level of theory increases from RHF to MP2 to CCSD(T)?
8.) Discuss your findings with respect to variations in basis set with a fixed level of theory, CCSD(T). Are there
trends in the percent errors as the basis set increases from aug-cc-pVDZ to aug-cc-pVTZ to aug-cc-pVQZ?
12
PART B: POTENTIAL CURVE, SPECTROSCOPIC CONSTANTS, AND DISSOCIATION ENERGY
OF DIATOMIC HYDROGEN (20 points)
In this part of the project, you will explore the ground potential energy curve of the diatomic molecule H2, extract
vibration-rotation spectroscopic constants, and determine the dissociation energy.
The chemical reaction for dissociation of a homonuclear diatomic molecule X2 into its constituent atoms in the gas
phase is given by
X2 (g)
→
2 X (g) .
(4)
The dissociation energy De is defined as the energy of the products minus the reactants in the reaction above,
De = 2E ( X ) − E ( X 2 ),
(5)
€
where E is the total electronic energy of the atom or molecule.
A typical potential energy curve of a diatomic molecule X2 which exhibits the bond dissociation into constituent X
atoms is shown in Figure 6. At low energies, the curve is approximately parabolic (i.e., the harmonic oscillator
approximation is valid), but as the energy increases the anharmonicity becomes more important.
Figure 6. Typical potential curve for a diatomic molecule, where R − Re is the bond displacement.
Beyond the harmonic oscillator, rigid rotor approximation, the quantized energy levels EvJ of a vibrating, rotating
diatomic molecule may be fit to a double polynomial expansion in the vibrational and rotational quantum numbers,
v and J. This is known as a Dunham expansion, which has the following form in wavenumbers,
ω vJ = ! 1$
! 1$
! 1 $2
2
EvJ
= ν! 0 # v + & + Be J ( J +1) − α e # v + & J ( J +1) − ν! 0 xe # v + & − D () J ( J +1)*+ + … .
" 2%
" 2%
" 2%
hc
(6)
Here, ω vJ is the vibration-rotation energy expressed in units of wavenumbers, ν! 0 is the harmonic vibrational
frequency in wavenumbers, Be is the rotational constant, α e is the vibration-rotation coupling constant, xe is the
anharmonicity constant, and D is the centrifugal distortion constant.
13
The spectroscopic constants in the Dunham expansion of the vibration-rotation energy may be determined from a
high-accuracy potential energy curve. Such a procedure can be found in Halpern, A. M. J. Chem. Ed. 2010, 87, 174179, for example. In this part of the project, you will obtain a potential energy curve for H2 from high-level
quantum mechanical calculations and fit it to a 6th order polynomial in the bond displacement in order to extract the
vibration-rotation spectroscopic constants as well as the dissociation energy. The level of theory you will employ is
QCISD (quadratically-convergent Configuration Interaction with double excitations); because H2 has only two
electrons, in this case QICSD is equivalent to full CI.
Procedure
1.) The first step in the analysis is obtain the equilibrium bond length and total electronic energy of H2 from a
geometry optimization calculation. Using Avogadro, build the H2 molecule. To begin, start Avogadro and click
on the Drawing Tool,
. Set the Element to Hydrogen, Bond Order to Single, and make sure that "Adjust
Hydrogens" is turned OFF. Click with your mouse in the drawing window to place a single hydrogen atom.
Then click again to place another hydrogen atom nearby (the H2 bond length is less than 1 Å). Use the mouse to
draw a bond between the hydrogen atoms in order to complete the construction of the H2 molecule.
Next, select "Extensions → Gaussian", and set up the calculation as follows:
Calculation: Geometry Optimization
Theory: RHF
Charge: 0
Format: Z-matrix
Processors: 4
Basis: 6-31G(d)
Multiplicity: 1
You should see on the second line in the white input box the following text:
#n RHF/6-31G(d) Opt
Edit the line so that it has the following form (you are replacing the level of theory and basis set and adding
'=Z-matrix' to the Opt keyword):
#n QCISD/aug-cc-pVQZ Opt=Z-matrix
In addition, you should add a line at the beginning of the white input box that has the form:
%mem=1GB
Once you have made all the changes in the settings, click "Compute" to submit the calculation. Once the
calculation is complete, close Avogadro and open the H2 results file with the 'gedit 'command. Scroll to the
bottom of the file, and look for a block of text similar to the following:
1\1\GINC-FRODO\FOpt\RQCISD-FC\Aug-CC-pVTZ\H2\STANDARD\26-Mar-2015\1\\#n QCISD
aug-cc-pVTZ OPT=Z-matrix\\H2\\0,1\H\H,1,R\\R=0.74309563\\Version=AS64LG09RevD.01\State=1-SGG\HF=-1.1330032\MP2=-1.1650011\MP3=-1.1705591\MP4D=1.1722869\MP4DQ=-1.1719563\MP4SDQ=-1.1720217\QCISD=-1.1726356\RMSD=9.864e09\RMSF=4.652e-05\Dipole=0.,0.,0.\PG=D*H [C*(H1.H1)]\\@
The total electronic energy is listed in bold above (QCISD=xxxxx) and is given in hartrees. Record this value
(your results will be slightly different than those listed above).
14
To obtain the equilibrium bond distance, search in the Gaussian results file for the string 'Stationary point
found'. The equilibrium value of the H2 bond distance is reported in angstroms a few lines below this string,
and just below that with more significant digits as part of the Z-matrix listing.
2.) The next step is to generate a portion of the potential energy curve of H2, which you will use to obtain
spectroscopic constants. Using Avogadro, build another copy of the H2 molecule. Select "Extensions →
Gaussian", and set up the calculation as follows:
Calculation: Single Point Energy
Theory: RHF
Charge: 0
Format: Z-matrix
Processors: 4
Basis: 6-31G(d)
Multiplicity: 1
You should see on the second line in the white input box the following text:
#n RHF/6-31G(d) Opt
Edit the line so that it has the following form (you are replacing the level of theory and basis set and changing
the calculation type to SCAN):
#n QCISD/aug-cc-pVQZ SCAN
In addition, you should add a line at the beginning of the white input box that has the form:
%mem=1GB
Finally, the last line in the white input box should look something like the following:
B1 0.75021
The numerical value may be different because it corresponds to the H-H bond length in the initial structure that
you built. Edit the line so that it has the following form:
B1 0.50 14 0.05
This sets up a scan of the potential energy curve of H2 which starts with a bond distance of 0.50 Å and computes
the total electronic energy of H2 at 14 additional distances using an increment of 0.05 Å.
Once you have made all the changes in the settings, click "Compute" to submit the calculation. When the
calculation completes, close Avogadro and open the results file using 'gedit'. Search for the string 'Summary of
the potential surface scan'. You should see a table of results that looks similar to those shown in Table 1. The
second column gives the H-H bond distance in Å and note that the desired QCISD energies are listed in the last
column of the results.
15
Table 1. Example of results from potential surface scan of H2. Note that these results were obtained at the
QCISD/aug-cc-pVTZ level, so they will be slightly different than the results from the QCISD/aug-cc-pVQZ level.
Summary of the potential surface scan:
N
B1
SCF
MP2
---- --------- ----------- ----------1
0.5000
-1.06343
-1.09475
2
0.5500
-1.09548
-1.12689
3
0.6000
-1.11549
-1.14700
4
0.6500
-1.12691
-1.15856
5
0.7000
-1.13214
-1.16397
6
0.7500
-1.13289
-1.16492
7
0.8000
-1.13036
-1.16263
8
0.8500
-1.12545
-1.15801
9
0.9000
-1.11881
-1.15170
10
0.9500
-1.11093
-1.14420
11
1.0000
-1.10217
-1.13586
12
1.0500
-1.09281
-1.12697
13
1.1000
-1.08307
-1.11775
14
1.1500
-1.07310
-1.10834
15
1.2000
-1.06333
-1.10027
---- --------- ----------- -----------
MP3
-----------1.09974
-1.13195
-1.15217
-1.16384
-1.16939
-1.17050
-1.16840
-1.16399
-1.15791
-1.15067
-1.14261
-1.13403
-1.12514
-1.11610
-1.10788
-----------
MP4DQ
-----------1.10078
-1.13305
-1.15332
-1.16507
-1.17070
-1.17191
-1.16992
-1.16563
-1.15970
-1.15261
-1.14473
-1.13634
-1.12765
-1.11883
-1.11086
-----------
QCISD
-----------1.10111
-1.13343
-1.15376
-1.16558
-1.17129
-1.17260
-1.17074
-1.16659
-1.16082
-1.15392
-1.14627
-1.13813
-1.12973
-1.12124
-1.11365
-----------
You will need to copy the table from the Gaussian results file into a separate text file which will then be used to
import the results into Microsoft Excel. To do this, start Microsoft Word on the Macintosh computer and open
a new blank document. Copy the rows corresponding to the potential surface scan results from the Gaussian
file (like those above) into this blank document. [To copy from the results file, select the desired rows of data
with the mouse and use Control-C. Copy into Word using Command-V.] Save the Word document on the
Desktop as a Plain Text File using "Save As", and give the file a descriptive name with the file extension '.txt'
(the default settings for Text Encoding and Line Breaks are fine). You should then copy the plain text file to a
removable drive for later analysis.
3.) Next, to determine the dissociation energy of H2, you will need to calculate the energy of the H atom using
Avogadro and Gaussian. To do this, follow the procedure from Part A#1 to construct the H atom. Select
"Extensions → Gaussian", and make sure that the following settings are entered for the calculation:
Calculation: Single Point Energy
Theory: RHF
Charge: 0
Processors: 4
Basis: 6-31G(d)
Multiplicity: 2
You should see on the second line in the white input box the following text:
#n RHF/6-31G(d) SP
Edit the line so that it has the following form (you are replacing the level of theory and basis set):
#n QCISD/aug-cc-pVQZ SP
Submit the calculation, and when it completes close Avogadro and open the results file using 'gedit'. Locate and
record the total electronic energy of the H atom.
16
Results, Analysis, and Discussion
1.) From the total electronic energies of the H2 molecule and H atom from steps #1 and #3 of the procedure,
compute the bond dissociation energy of H2 at the QCISD/aug-cc-pVQZ level of theory using Equation (5).
Report your result in both atomic units and cm–1 [the conversion factor is 1 a.u. = 219474.6 cm–1].
2.) Compare your results for the equilibrium bond length and dissociation energy of H2 at the QCISD/aug-cc-pVQZ
level of theory with experimental literature values. Compute percent differences. Which agrees better with
experiment, the equilibrium bond length or the dissociation energy? You also should compare your QCISD
results to results from other lower levels of theory; some results may be found in Table 11-9 of Lowe and
Peterson. Discuss the degree of variation in the bond length and dissociation energy as the level of theory is
improved.
3.) Import your H2 potential energy curve data at the QCISD/aug-cc-pVQZ level from step #2 of the procedure into
Microsoft Excel. Include a table listing this data in your project.
Next, calculate bond displacements by subtracting the equilibrium bond distance from each value of the bond
distance; convert these bond displacement values into units of bohr. Produce a graph with the QCISD energy
(in hartrees) on the y-axis and bond displacement (in bohr) on the x-axis. Include this graph in your project.
Fit the data in the graph to a 6th-order polynomial using Microsoft Excel. Include a plot with this trendline in
your project and make sure that the equation and R2 value are displayed on the chart. It is also helpful to
modify the trendline format to display the parameters of the fit to about 5 or 6 decimal places.
The 6th-order polynomial representing the potential energy V of H2 has a functional form given by
2
3
4
V (R) = A0 + A1 ( R − Re ) + A2 ( R − Re ) + A3 ( R − Re ) + A4 ( R − Re ) 5
6
(7)
+ A5 ( R − Re ) + A6 ( R − Re ) .
In this equation, R is the bond distance, Re is the equilibrium bond distance (and thus R − Re is the bond
displacement), and A0 through A6 are the parameters of the fit. Note that since the polynomial expansion is
centered about the equilibrium bond distance (i.e., the minimum of the potential curve), the parameter A1
should be close to zero, though usually not exactly zero due to round-off errors and such.
The dissociation energy De for H2 may be determined from Equation (5) using the total electronic energy of the
H atom from step #3 of the procedure, and the total electronic energy of H2 as parameter A0 , which should
closely match the total electronic obtained in step #1 of the procedure. Compute and report the dissociation
energy De of H2. Note that the value of De obtained here should closely match the value obtained in #1,
though there may be small differences due to the fitting procedure.
The harmonic vibrational frequency ν 0 is given by
ν 0 = €
1
2π
! k $1/2
# & ,
"µ %
(8)
where k is the harmonic force constant and µ is the reduced mass of the diatomic molecule,
µ =
€
m1 m2
,
m1 + m2
(9)
17
where m1 and m2 are the masses of the two atoms (in this case, both are hydrogen atoms). We showed in class
that when the potential energy function V (R) is expanded in a polynomial series in terms of the bond
displacement, the harmonic force constant equals two times the parameter of the second-order term (i.e., the
term has the form 12 kx 2 , where x is the bond displacement). In this case, therefore, the force constant k is
defined as
k = 2A2 .
(10)
Use this relationship to obtain the force constant k from your fit, (note that A2 has units of hartree bohr–2) and
then use Equation (8) to determine the harmonic vibrational frequency ν 0 in s–1. Make sure to use a reduced
mass with at least five significant figures so that you get an accurate value for the harmonic frequency. Convert
into units of cm–1 and report the results for the harmonic frequency using the relation
ν0
,
c
ν˜ 0 =
(11)
where ν˜ 0 is the harmonic frequency in units of cm–1.
€
The rotational constant Be can be determined from the equilibrium bond length and reduced mass using the
€ relation
Be = h
.
8π cµ Re2
2
(12)
Compute the rotational constant Be for H2 and report its value in cm–1.
The centrifugal distortion constant D (not to be confused with the dissociation energy De ) is defined by the
equation,
D = 4Be3
.
ν 02
(13)
Using your calculated values of the harmonic vibrational frequency ν! 0 and rotational constant Be in cm–1,
determine and report the centrifugal distortion constant D in cm–1.
The vibration-rotation coupling constant α e may be determined from the following relation,
α e = −
6Be2 "
AR %
$1 + 3 e ' .
ν! 0 #
A2 &
Using Re in bohr, the previously-calculated values of the harmonic vibrational frequency ν! 0 and rotational
constant Be in cm–1, along with the parameters A2 and A3 from the 6th-order fit, calculate and report the
vibration-rotation coupling constant α e in cm–1.
(14)
18
Finally, the anharmonicity constant xe may be determined using the equation,
xe = 2
(
+
!
Be *
α ν! $
A R2
15 #1 + e 20 & − 12 4 e - .
8ν! 0 *)
A2 -,
6Be %
"
(15)
In this equation, if α e , Be , and ν! 0 are expressed in cm–1, with Re in bohr, then xe will be unitless. Using the
previously-calculated spectroscopic constants, along with the fit parameters A2 and A4 , calculate and report
the anharmonicity constant xe .
In summary, when you have completed this section, you should have a graph of the energy vs. bond
displacement for H2 along with a 6th-order polynomial fit of the data. In addition, you should have determined
the dissociation energy De along with the following five vibration-rotation spectroscopic parameters for H2: (1)
the harmonic vibrational frequency ν! 0 , (2) the rotational constant Be , (3) the centrifugal distortion constant D,
(4) the vibration-rotation coupling constant α e , and (5) the anharmonicity constant xe . These items all should
be included in your project. It might be helpful to show your work for any calculations in order to receive
partial credit in the case that something is in error.
4.) Obtain experimental literature values for each of the five vibration-rotation spectroscopic constants. Cite your
source(s). Compare your results to the literature, being quantitative in your comparisons. Can you suggest an
approach to obtain even better agreement between the calculations and experiment?
19
PART C: COMPLETE BASIS SET EXTRAPOLATION OF THE ENERGY AND PROPERTIES OF WATER
(15 points)
In this part of the project, you will explore quantum mechanical calculations on the water molecule, H2O. Because
water is a small molecule with only one non-hydrogen atom, very high-level calculations with large basis sets may
be performed in order to obtain highly accurate energies and properties.
The calculations will employ the CCSD(T) level of theory, which has been referred to as the "gold standard" of
quantum chemistry for molecular systems [Cramer, C. J. Essentials of Computational Chemistry, Wiley, New York,
2002, pp. 211-213]. This refers to the fact that when the basis set size is increased, energies and many properties
computed at the CCSD(T) level of theory appear to converge to asymptotic values, known as the complete basis set
(CBS) limit. And for many systems, the CBS extrapolation of the CCSD(T) results provides excellent agreement
with experimental results. Please refer to the many examples of CBS extrapolation found in Dunning, T. H. J. Phys.
Chem. A 2000, 104, 9062-9080.
There are a variety of functional forms used to obtain the asymptotic limit from the CBS extrapolation. One of the
most comment is the X–3 fit, in which the property of interest P is calculated using correlation-consistent basis sets,
aug-cc-pVXZ, where X=D, T, Q, 5, etc. The values of the property P for each basis set level, P ( X ) , is then fit to
the following functional form,
P ( X ) = PCBS + A
.
X3
(16)
In this equation, X=2, 3, 4, 5. Also, PCBS is the value of property P in the limit of a complete basis set (i.e., as
X → ∞ ) and A is a fitting parameter. A typical CBS extrapolation is shown in Figure 7 for the determination of the
dissociation energy of N2 at the MP2, MP3, and MP4 levels of theory with cc-pVXZ (X=D, T, Q, 5) basis sets.
Figure 7. CBS extrapolation of the dissociation energy of N2 at various levels of theory [from
Dunning, T. H. J. Phys. Chem. A 2000, 104, 9062-9080].
In this part of the project, you will perform CBS extrapolations of the energy, geometry, and dipole moment of the
water molecule in order to investigate their asympototic convergence and obtain these properties to high accuracy.
The agreement between the CBS-extrapolated results and experimental values will be explored.
20
Procedure
1.) The first thing you need to do is build the H2O molecule using Avogadro. To construct an initial structure for
the H2O molecule, start Avogadro and click on the Drawing Tool,
. Set the Element to Oxygen, Bond
Order to Single, and make sure that "Adjust Hydrogens" is turned ON. Click with your mouse in the drawing
window to place a single oxygen atom. This should actually place an oxygen atom along with two attached
hydrogens in the drawing window (i.e., you have built the H2O molecule).
Next, perform your first quantum mechanical calculation of the H2O molecule using the Gaussian software
package. Select "Extensions → Gaussian". Use the following settings:
Calculation: Geometry Optimization
Theory: RHF
Charge: 0
Processors: 4
Basis: 6-31G(d)
Multiplicity: 1
You should see on the second line in the white input box the following text:
#n RHF/6-31G(d) Opt
Edit the line so that it has the following form (you are replacing the level of theory and basis set and adding
'=Z-matrix' to the Opt keyword):
#n QCISD/aug-cc-pVQZ Opt=Z-matrix
In addition, you should add the following line to the beginning of the white input box:
%mem=1GB
Once you have made all the changes in the settings, click "Compute" to submit the calculation. When the
calculation is complete, close Avogadro and open the H2 results file with the 'gedit 'command. You will need to
extract a variety of information from the results file, including the total electronic energy, CPU time, number of
basis functions, and equilibrium geometry.
To obtain the total electronic energy, scroll to the bottom of the file, and look for a block of text similar to the
following:
1\1\GINC-FRODO\POpt\RCCSD(T)-FC\Aug-CC-pVDZ\H2O1\STANDARD\26-Mar-2015\1\\#n
CCSD(T) aug-cc-pVDZ OPT=Z-matrix\\H2O\\0,1\H\O,1,B1\H,2,B2,1,A1\
\B1=0.96669923\B2=0.96670035\A1=103.87253635\\Version=AS64L-G09RevD.01\State=1A'\HF=-76.040692\MP2=-76.2609086\MP3=-76.265512\MP4D=-76.2695628\MP4DQ=76.2670567\MP4SDQ=-76.2688246\CCSD=-76.2686226\CCSD(T)=-76.2739035\RMSD=2.317e09\RMSF=1.559e-04\PG=CS [SG(H2O1)]\\@
The total electronic energy is listed in bold above (CCSD(T)=xxxxx) and is given in hartrees. Record this
value. You should also record the value of the energy computed at the Hartree-Fock level (HF=xxxxx) for later
use in exploring the electron correlation effects.
To obtain the dipole moment, look in the results file just above the block of text with the energy. You should
see some lines listing the dipole moment, such as shown below:
Dipole moment (field-independent basis, Debye):
X=
0.0000
Y=
2.0159
Z=
0.0000
Tot=
2.0159
The value you want is the total, listed at the end of the line (in this case, 2.0159); the units are Debye.
21
Before you close the results file, you also should record the CPU time to the nearest minute, which is listed at
the very end of the file after the random quote. And finally, you should record the number of basis functions,
which may be found by searching for the string "Standard basis". You should find a block of text similar to the
following:
Standard basis: Aug-CC-pVDZ (5D, 7F)
There are
20 symmetry adapted cartesian basis
There are
4 symmetry adapted cartesian basis
There are
7 symmetry adapted cartesian basis
There are
12 symmetry adapted cartesian basis
There are
18 symmetry adapted basis functions
There are
4 symmetry adapted basis functions
There are
7 symmetry adapted basis functions
There are
12 symmetry adapted basis functions
41 basis functions, 65 primitive gaussians,
functions of A1 symmetry.
functions of A2 symmetry.
functions of B1 symmetry.
functions of B2 symmetry.
of A1 symmetry.
of A2 symmetry.
of B1 symmetry.
of B2 symmetry.
43 cartesian basis functions
The total number of basis functions is listed in bold. Record this value.
To obtain the equilibrium bond distances and bond angle, seach in the Gaussian results file for the string
'CONVERGENCE CRITERIA'. The equilibrium values of the H2O bond distances (in angstroms) and bond
angle (in degrees) are reported a few lines below this string under 'Optimized Parameters', and just below that
with more significant digits as part of the Z-matrix listing. Record the bond distances to 5 decimal places (they
should be equivalent to that number of digits) and the bond angle to 4 decimal places.
2.) Next, you should repeat the calculation of water for the TZ, QZ, and 5Z basis sets. Follow the same procedure
as in step #1, except perform three more calculations in which you change the basis set to the following:
a) aug-cc-pVTZ
b) aug-cc-pVQZ
c) aug-cc-pV5Z .
All the other settings should be as in #1. The aug-cc-pVTZ calculation may be run as usual by clicking
"Compute" from the Avogadro window to submit the calculation. The aug-cc-pVQZ and aug-cc-pV5Z
calculations are more lengthy and hence you should submit these using a different method so that you can have
the calculations run while you are off doing other things.
For the QZ and 5Z calculations, from the Avogadro setup window for Gaussian, click "Generate". You will be
prompted to enter a filename. Make sure to record the name of the file for future reference; it should end with a
'.com' extension.
To submit the calculation, close Avogadro, and in the terminal window type the following command:
g09
your-h2o-file.com
your-h2o-file.log
&
Here, 'your-h2o-file.com' corresponds to the name of the file that you saved with the "Generate" command. The
command above submits your Gaussian calculation. Run one calculation at a time. Do not try to run both the
QZ and 5Z calculations simultaneously; it will bog down the computer significantly. The QZ calculation
should take about 30 minutes or so of real time, while the 5Z calculation should take about 6 hours of real time
(note that since you are using 4 processors for the calculations, the amount of CPU time is roughly 4× the real
time). Once you have submitted a calculation, you may log off and return later to get the results. Make sure
that you log onto the same Linux computer where you submitted the calculation (i.e., if you submitted the QZ
calculation on frodo, make sure to go back to frodo to get the results).
22
Once the TZ, QZ, and 5Z calculations are finished, follow the instructions in step #1 of the procedure to record
the HF and CCSD(T) total electronic energies, dipole moment, CPU time, number of basis functions, and
equilibrium bond distances and angle for each calculation.
Results, Analysis, and Discussion
1.) Include in your project a table listing the basis set, number of basis functions, CPU time, HF energy, CCSD(T)
energy, dipole moment, bond distance, and bond angle for H2O at the CCSD(T) level with the aug-cc-pVXZ
(X=D, T, Q, and 5) basis sets.
2.) For the CCSD(T) energy, dipole moment, bond distance, and bond angle, make a graph with the property on the
y-axis and the ordinal number of the basis set (X=2, 3, 4, 5) on the x-axis. You should have four graphs total;
turn these in with your project. Discuss whether or not each of the properties appears to be approaching an
asymptotic value.
3.) To fit the properties to the asymptotic form given in Equation (16), for those properties studied in #2 above
create a graph with the property on the y-axis and 1/X3 on the x-axis. Only do this for those properties that
appear to go to an asympotote. You should have a minimim of two graphs total for this part; turn these in with
your project. The DZ results often deviate somewhat from the trend and are therefore usually left out of the fit.
Therefore, you should fit only the X=3,4, and 5 data points to a linear trendline (this is referred to as a TQ5 fit);
make sure to include the equation of the line and the R2 value on the graph. Also, even though you are not
including the X=2 data points in the fit of the trendline, make sure to show them on the graphs; you may include
the X=2 data point as a separate series on each graph if necessary.
4.) From the trendlines for each graph in #3, tabulate and report the values of PCBS and A for each of the properties
with trendlines. While there is no experimental value of the total electronic energy at the CCSD(T) level, you
should be able to find experimental literature values of the gas phase dipole moment, bond distances, and bond
angle of H2O (cite your sources).
Prepare a table of results listing the values of the properties (dipole moment, bond distance, and bond angle) for
each basis set and at the CBS limit if available. Compare your results at each basis set level and at the CBS
limit with the literature value; use percent differences.
Discuss your findings. Are there properties which show little or no improvement as the basis set increases? Are
there others that show improvement? For those properties for which you have a CBS limit, discuss the
agreement between the CBS results and experiment. Do you think that the CCSD(T)/CBS method holds up as a
"gold standard" in these cases?
5.) The difference between the exact electronic energy and the Hartree-Fock energy is defined to be the correlation
energy, Ecorr . The correlation energy may be estimated by assuming that the CCSD(T) energy is close to the
exact energy (i.e., we assume that the CCSD(T) level recovers the vast majority of the electron correlation
effects). Thus, we can calculate the correlation energy using the equation
Ecorr ≈ E (CCSD(T )) − E ( HF ) .
(17)
Use the HF and CCSD(T) total electronic energies from the aug-cc-pVXZ (X=D, T, Q, and 5) calculations to
calculate a value of Ecorr for the DZ, TZ, QZ, and 5Z basis sets. Report these values.
23
In addition, you already have the CBS energy at the CCSD(T) level from step #4. To determine the CBS
energy at the HF level, use the HF energies from the aug-cc-pVXZ (X=D, T, Q, and 5) calculations and the
same procedure for fitting to the X–3 form that you employed in step #3 to obtain ECBS ( HF ) . Use this to
calculate the correlation energy, Ecorr , from the CBS-extrapolated CCSD(T) and HF energies. Report this
value.
Compare your results at the DZ, TZ, QZ, 5Z, and CBS levels for Ecorr for H2O with the "exact" value reported
in Table 11-4 of Lowe and Peterson. Discuss the agreement. Compute the percentage of total correlation
energy recovered in each case by taking the ratio of the calculated correlation energy to the exact correlation
energy (from Lowe and Peterson). Based on the results presented in Figure 3 of Dunning (2000), assuming
most of the correlation energy of H2O comes from the oxygen atom, how well do your %correlation values
match with the literature?
6.) From your CPU times and basis set sizes, make a graph with the natural log of CPU time (in minutes) on the yaxis and natural log of the basis set size on the x-axis. Fit the data to a linear trendline and display the equation
and R2 value on the chart; include this graph in your project.
The graph allows you to determine the approximate scaling of the CCSD(T) method with basis set size. In
general, the CPU time of a method may be given by the following relation,
CPU time = cK S ,
(18)
where K is the number of basis functions, c is a proportionality constant, and S is the scaling of the method.
This relationship may be converted into a linear one by taking the natural log of both sides to yield,
ln (CPU time) = ln ( c ) + S ln ( K ) .
(19)
Here we see that a graph like the one you created should be linear with a slope corresponding to the scaling
factor S.
For the CCSD(T) method in the Gaussian software package, what do you find to be the scaling? Report S.
Based on your computed scaling factor S, estimate (and report) how much CPU time (in minutes, then convert
to hours and days) it would take to perform a CCSD(T)/aug-cc-pV5Z calculation on ethanol, CH3CH2OH.
Using the aug-cc-pV5Z basis set, the basis set size for ethanol is 861 functions.