Quantum Chemistry tutorial


Exploring homolytic bond cleavage of a hydrogen molecule with Gaussian09


Introduction

h2 In the previous lecture we've learned about the theoretical basis of electronic structure methods, and in particular about the Hartree-Fock method. Now it is time to jump into the electronic structure business and perform a couple of calculations by yourself. The package we are going to use is called
Gaussian. It is among the modern electronic structure codes available and provides an easy-to-use interface called Gaussview that allows for a user-friendly access to quantum chemistry. Additionaly the Gaussian program provides many post-Hartree-Fock methods.

To learn how to use the program by yourself or maybe extend your knowledge beyond the scope of this tutorial one can check out the links on the following site.

To start a calculation we basically need four things:

The files needed for this practical can be downloaded as an archive here and unpacked by typing

tar xzvf h2.tar.gz
However, you don't need these, as all files required can be created with the information below.

The H-H molecule

Our calculations will be on the simple H2 molecule, which is the smallest molecule that exhibits all fundamental aspects of chemical bonding (why?).

Single point energy calculation

A simple input for the electronic structure calculation of H2 is given below:
%chk=h2_SP
#P RHF 6-31G 

h2 molecule

0 1
h 
h 1 r

r 1.0

The first line starting with %chk specifies that we want to keep the checkpoint file. You are free to choose a file name; I used h2_SP (standing for single Point). Orbitals are stored on this file and can be visualized with Gaussview. The #P tells Gaussian to print all information during the computation. Then, the method is specified, which is Restricted Hartree Fock (RHF) and the basis (6-31G) is selected, using the so-called Pople notation. G stands for "Gaussian", so a GTO basis with 6 gaussian functions for the inner (non-valence) electrons and 4 (3+1) for the valence electrons is used. Then comes a whileline followed by a title string and another whiteline. The next two integers represent the total charge (0 in this case), and spin multiplicity (i.e. the value of 2*S+1, where S is the eigenvalue of the S*S operator: S(S+1)). Since we have a closed shell system, we have S(S+1)=0; thus S=0 and 2S+1 = 1, a Singlet. For the geometry we use a Z-matrix. In the case of a diatomic molecule, one needs to specify only the interatomic distance r, which we set at 1.0 Angstrom.

Now we can do the electronic structure calculation by typing

g09 h2.com 
During the calculation Gaussian writes an output file that communicates the most important infos to the user. You can extract the following information: Additionaly you could extract the individual HF energies after each SCF cycle and plot them against the number of cycles. I assume you know how at least one program to visualize data

Molecular orbitals

In the next step we want to see the bonding and antibonding molecular orbitals (MO) of the H2 molecule. For this we have to start the Gaussview program by typing
gview.exe 
In the main window, select File to open the h2_SP.chk checkpoint file. The click the molecular orbital button

to open the orbital interface. There, select the Visualize tab and click Update to render the orbitals.

To see the orbtials, click on the small square in the energy diagram:

Dissociation into H atoms

So far, we have only performed so-called single-point calculations, i.e. computed the HF energy for a fixed nuclear configuration. Gaussian is quite versatile when it comes to automatically perform scans along selected coordinates. In our case of the H2, we will now determine the energy profile (or the potential energy surface (PES)) along the H-H distance. The input file now looks like this:
%chk=h2_rhf_scan
#P RHF 6-31G SCAN

h2 molecule

0 1
h
h 1 r

r 0.5 S 55 0.1

The Scan keyword specificies that we want to scan one or more coordinates. In our case, we scan along the intermolecular distance, starting from 0.5 angstrom up to 6.0 Angstom in 55 steps of 0.1 Angstrom each. If you think this takes too long to compute, you can reduce the number of steps, and increase the stepsize accordingly. The program will loop over all specified distances and perform for each of those an independent SCF calculation. Finally, a table of the respective HF energies will be written to the output (.log) file. Gaussview can be used to create plots, but you're free to use a program of choice. In Gaussview, open the .log file, select in the main window the Results menu and choose Scan. This will bring up a new window showing the electronic energy as a function of the H-H distance.

Look at the curve. What is the energy one would expect for the dissociation limit (i.e. r-> oo) of H-H?

As will be discussed in the lectures, at the dissociation limit, the restricted Hartree Fock solution contains so-called ionic contributions, in which both electrons are localized on one of the atoms. This leads to a too high dissociation energy (by how much?). Clearly these contribution should decrease with increasing intermolecular distance r.

The Unrestricted Hartree Fock method, in which the paired alpha and beta spins are not forced to occuppy the same orbital, overcomes this problem (but wavefuntion is not eigenfunction of the S*S operator). Use the UHF method and repeat the dissociation calculation.

The input file should look like

%chk=h2_uhf_scan
#P UHF 6-31G SCAN guess=(mix,always)

h2 molecule

0 1
h
h 1 r

r 0.5 S 55 0.1

Please observe that a broken symmetry UHF singlet wavefunction is only obtained using the guess=(mix,always) keyword. If this is not used, the initial guess is chosen such that the SCF calculation converges to the RHF wavefunction even with the UHF Ansatz.

Beyond Restricted Hartree Fock

As pointed out in the lecture, the HF method gives molecular interaction energies that are only approximations to the exact quantum many-body energy. In the following we will study the effect of improving the RHF approximation. This will be achieved via perturbation theory, the so-called Moeller-Plesset (MPx, where x gives the order in the perturbation expansion), multi-configurational methods (truncated CI), coupled cluster and variational multi-determinant methods, like e.g. CASSCF.

We will create an input file for each of these methods and repeat the scans.

For MP2

%chk=h2_mp2_scan
#P MP2 6-31G SCAN

h2 molecule

0 1
h
h 1 r

r 0.5 S 55 0.1

For MP4
%chk=h2_mp4_scan 
#P mp4 6-31G SCAN

h2 molecule

0 1
h
h 1 r

r 0.5 S 55 0.1

Because there are only two electrons, truncating the configuration interaction expansion after doubles, is same as full CI for this system.
%chk=h2_cisd_scan
#P CISD=(MaxCyc=100) 6-31G SCAN

h2 molecule

0 1
h
h 1 r

r 0.5 S 55 0.1

In the lecture it was mentioned that truncating the CI expansion at the singles level (i.e. only consider single electron excitation from the occuppied into the unoccuppied orbitals) does not correlate with the ground state. NOw you can check this yourself:
%chk=h2_cis_scan
#P CIS=(MaxCyc=100) 6-31G SCAN

h2 molecule

0 1
h 
h 1 r

r 0.5 S 55 0.1

For CASSCF. Note that since we include all electron (2) and all orbitals (4), it is essentially full CI, apart form fact that the MO coefficients are optimized. See if this affects the outcome.
%chk=h2_cas_scan
#P CASSCF(2,4) 6-31G SCAN

h2 molecule

0 1
h
h 1 r

r 0.5 S 55 0.1

However, the idea of the CASSCF method is to limit the CI expansion to only those orbital that matter (or that we believe matter). Try reducing the number of orbitals to include only the 1-sigma bonding and anti-bonding orbitals. Finally, we try the dissociation with the Coupled-Clusters-singles-doubles (CCSD) method.
%chk=h2_CCSD_scan
#P CCSD 6-31G SCAN

h2 molecule

0 1
h
h 1 r

r 0.5 S 55 0.1

After completing these calculations you can compare the dissociation curves. In our example the CISD method is same as full CI, and hence be considered as the most accurate. Compare your result with respect to the obtained equilibrium bonding distance (How large is the error for this quantity?) and the dissociation limit for each method. Hereby the CISD results should serve as the reference energy curve.

Compare the results of the calculations with second and fourth order perturbation (MP2 -> MP4)

What is the error in the CI result compared to the analytical results of -1 Hartree. How can we improve the CI dissociation energy?

The final results are summarized in the figure below:

This ends the tutorial. However, the H2 molecule is too simple to illustrate the performance of the ab initio quantum chemistry tool box. Therefore, if you have time, you may want to try out a more complicated di-atomic molecule, such as F2, N2, CO, O2. In the latter cases, keep in mind that the spin state is triplet (why).