Correct treatment of bond dissociation

Correct treatment of bond dissociation
Theoretical chemistry profile course – computer exercise
Purpose
The purpose of the exercise is to:
•
•
•
test the limitations of common quantum chemistry methods
get familiar with the MOLCAS software package
get experience with complete active space (CAS) calculations
Introduction – bond dissociation
Most quantum chemistry methods are relatively good at describing closed shell systems
and equilibrium structures. However, the description of bond dissociation is much more
challenging. This can be important e.g., in photochemistry when excitation of a molecule
leads to photodissociation, or in chemical reactions where the reaction goes through a
transition state where the bonds are elongated relative to the equilibrium structures.
In the diagnostic test you calculated the potential energy profile for the H2
molecule using restricted Hartree-Fock (RHF). The calculation did not converge to the
correct dissociation limit and the energetic error was large. In this exercise you will study
the same reaction using a number of different methods, which should give you some
insight into how these methods describe bond dissocation. It is recommended that you do
the calculations using the MOLCAS program, which will be used in some of the coming
exercises.
General procedure
Calculate the energy of H2 as a function of distance. Use the energy of two isolated H
atoms as a reference (zero of energy). Note that the energy even of the single-electron H
atom varies with the choice of functional in the DFT calculations.
1a. Repeat the calculations of singet and triplet state using HF. however use the cc-pvtz
basis set instead of the aug-cc-pvdz basis set you used in the diagnostic test.
1b. Calculate the dissociation of singlet H2 in UHF. Do you get the result that you
expect?
1c. Check the expectation value of the S2 operator S(S+1). What value should you get for
a real singlet? What is the value in the UHF calculation? What does this value say about
the wavefunction?
2. Perform the calculations using the DFT functional BLYP, both restricted and
unrestricted. What do the results say about the possibility to describe static correlation
with DFT functionals like BLYP?
3. Perform the calculations using at least one other method of your choice, both restricted
and unrestricted.
4. Perform the calculations using the CASSCF method. For this method you need to
select the number of active electrons that are correlated, and the number of active orbitals
(both occupied and virtual).
6. Select a suitable method to calculate the dissociation energy curve of the O2 molecule.
Useful MOLCAS inputs
For this lab, we’ll be using mostly MOLCAS. The input structure of MOLCAS is very
different from the Gaussian one. The input calls different modules (starting with &), each
of them having their respective input in the form " keyword = value ".
Single HF calculation
&GATEWAY
Coord
2
H 0.0 0.0 0.0
H 0.75 0.0 0.0
Basis Set=cc-pVTZ
Group=Nosym
&SEWARD
&SCF
Scan the H-H distance
To scan a bond distance the MOLCAS input uses loops:
>>> Foreach DISTR in (5 .. 50)
>>> Eval DIST = $DISTR/10
&GATEWAY
Coord
2
H 0.0 0.0 0.0
H $DISTR 0.0 0.0
Basis Set=cc-pVTZ
Group=Nosym
&SEWARD
&SCF
>>>End Do
This gives a scan of the H-H bond from 0.5 to 5.0 Å starting with steps of 0.1 Å.
Create symmetry breaking
Convergence of the SCF does not guarantee that the energy is a global minimum, in can
either be in a local minimum or a saddle point. You can artificially create symmetry
breaking (of amplitude 0.2) to leave the saddle points:
&GATEWAY
Coord
2
H 0.0 0.0 0.0
H 0.75 0.0 0.0
Basis Set=cc-pVTZ
Group=Nosym
&SEWARD
&SCF
UHF
scramble=0.2
CASSCF calculation
The main strength of MOLCAS is the multi-reference calculation abilities. A CASSCF
calculation normally requires two steps: the active space design and the actual
calculation. We will use this strategy in later labs. For H2, the active space is trivial (the σ
and σ*) and stable, so one can run directly:
&GATEWAY
Coord
2
H 0.0 0.0 0.0
H 0.75 0.0 0.0
Basis Set=cc-pVTZ
Group=Nosym
&SEWARD
&RASSCF
Inactive = 0
ras2 = 2
This means that there are no inactive (occupied not in the active space) orbital but 2
orbitals in the ras2 space (that is the CAS active space).
Running MOLCAS
To run MOLCAS, one simply runs:
molcas -f somename.input
The output file will then be created with name "somename.log", and plenty of other
important files will be generated (orbitals, optimized geometries, …).
Some of this files are meant to be visualized, for example the .molden files are visualized
with molden.
The MOLCAS output is often very readable (for the most common programs), but it is
still practical to simply grep the energies. Most common programs print the energy after
"::".
grep "::" somename.log
Report
Plot the dissociation energy curves for all methods. Answer the questions from the
instructions and add any additional observations that you find relevant.
Appendix: Testing wavefunction stability in Gaussian
In Gaussian, you would checking the stability of your SCF solution using the keyword
stable:
%chk=filename.chk
#p uhf/cc-pvtz scf=tight stable=opt
Comment line: H-H bond stretch
0 1
H
H 1 r1
r1 5.5
The keyword stable=opt tells the program to check the stability of the wavefunction,
and in case it is not stable, to optimize it to a new stable solution. The stability check
should ideally be performed for all calculations where there could be any doubt as to the
nature of the wavefunction. The inclusion of the command %chk=filename.chk tells
the program to save important information about e.g., the converged wavefunction/
density in the file filename.chk (use the same name as for the input file, but with the .chk
extension).
When continuing a calculation, it is recommended that you use the converged
wavefunction of from previous run as a starting point. You do that by first copying the
checkpoint file
>cp filename.chk filename2.chk
and then calling for the new file with guess=check. If you use guess=check the
program expects to find the file listed in the %chk line at the beginning of the
calculation. It will then overwrite the information in that file with the results from the
present calculation. Example input file:
%chk=filename2.chk
#p uhf/cc-pvtz scan scf=tight guess=check
Comment line: H-H bond stretch
0 1
H
H 1 r1
r1 5.5 50 -0.1
The calculation will scan from 5.5 Å towards lower energy, using the stable UHF
wavefunction as an initial guess.

Download Report