# Introduction to Electronic Structure Calculations using MOLPRO

#### Dr. Udo W. Schmitt, Max-Planck Institut fuer Biophysikalische Chemie,Goettingen

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 MOLPRO, which has been developed

by Peter Knowles and Hans-Juergen Werner over the last 15 years. It is among

the modern electronic structure codes available and provides an easy-to-use scripting language that allows for

a user-friendly access to quantum chemistry. Additionaly MOLPRO provides for the largest set of implemented 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 "Quick Start" or "User's Manual" on the following site

To start a calculation we basically need four things:

- the structure of the molecule of interest in (cartesian) coordinates
- the overall charge of our system
- the method that will be used (Hartree-Fock (HF),
Density Functional Theory (DFT), Moeller-Plesset perturbation theory (MPx) etc.)

- a basis set.

The files needed for this practical can be downloaded as an archive here

and unpacked by typing

tar xzvf prakt1.tar.gz

### The H-H molecule

Our first calculation will be the simple H2 molecule, which is the smallest

molecule that exhibits all fundamental aspects of chemical bonding (why?).

A simple input for the electronic structure calculation of H2 is given below:

basis=6-31G

r=1.0 ang

geometry={angstrom;

h1;

h2,h1,r}

hf

First, the basis (6-31G) using the so-called Pople notation is specified. 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, the H-H distance and the so-called internal coordinates (bonds,angles, torsions) are specified.Finally, the key word "hf" stands for "do a Hartree-Fock calculation"

Now we start out first electronic structure calculation by typing

molpro h2.com

During the calculation MOLRPO writes an output file that communicates the most

important infos to the user. Please extract the following information:

- number of SCF cycles needed for convergence
- the HF energy (what units are being used?)
- the dipole moment

Additionaly you could extract the individual HF energies after each SCF cycle and plot them

against the number of cycles.

In the next step we want to perform the visualization of the bonding and antibonding molecular

orbitals (MO) of the H2 molecule. For this we have to redo are electronic structure calculation:

molpro h2_orbitals.com

The same calculation as in the previous step will be performed, but additionaly the MOs

and the electron density are written to external files. To view the MOs, we will have to use VMD, a molecular visualization program for displaying, animating, and analyzing large biomolecular systems using 3-D graphics and built-in scripting.

Typing

vmd -e orbitals.vmd

will load an existing script that loads the relevant MOs.

In the VMD main window you can activate(deactivate) the two MO and the electron density while clicking on the D entry for each MO/density entry in the main window. Which orbital is the bonding, which one is the anti-bonding one?

So far, we have only performed so-called single-point calculations, i.e. computed the HF energy

for a fixed nuclear configuration. MOLPRO 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:

r=[0.3 0.4 0.6 0.8 1.0 1.2 1.4 1.6 2.0 4.0 6.0] ang

basis=6-31G

i=1

geometry={

h1;

h2,h1,r(i)}

do i=1,#r !loop over range of bond distances

hf

e1(i)=energy

enddo

{table,r,e1

title,HF PES for H-H

plot,file='h2_hf_pes.plot'

}

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 a file called *h2_hf_pes.plot*

Now we going to use the plot program "xmgrace" that comes with most Linux distributions

to look at the PES. Typing

xmgrace h2_hf_pes.plot

will reveal the electronic HF energy as a function of the H-H distance.

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

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 HF approximation. This will be achieved via perturbation theory, the so-called

Moeller-Plesset (MPx, where x gives the order in the perturbation expansion) and variational

multi-determinant methods, like e.g. configuration interaction (CI). In molpro this can simply be

achieved by adding the respective keywords:

r=[0.3 0.4 0.6 0.8 1.0 1.2 1.4 1.6 2.0 3.0 4.0 5.0 6.0] ang

basis=6-31G

i=1

geometry={

h1;

h2,h1,r(i)}

do i=1,#r !loop over range of bond distances

hf ;e_hf(i)=energy

mp2,NOCHECK ;e_mp2(i)=energy

casscf ;e_cas(i)=energy

mrci ;e_ci(i)=energy

enddo

{table,r,e_hf,e_mp2,e_cas,e_ci

title,HF PES for H-H

plot,file='h2_methods_pes.plot'

}

After a sucessful finish of the MOLPRO run you can have a look at the four PESs again using "xmgrace".

In our example the CI method (mrci) can be considered as the most accurate. CASSCF stands for "complete active space SCF", which can be regarded as a truncated CI expansion. 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 CI results should serve as the reference energy curve. Repeat the calculation with changing the

order in the perturbation expansion from 2 to 4 (mp2 -> mp4) and compare the two results. 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 water molecule

Now we going to move to a slightly more complex molecule having three nuclei and 10 electrons: the water monomer. For the structural information MOLPRO allows for the use of

the convient xyz in cartesian coordinates format. An example xyz-File is shown here:

3

remark line

O 0.00000 0.00000 0.00000

H 0.00000 0.00000 0.94000

H 0.00000 0.99030 -0.23192

With this structure, which does not represent the well-known equilibrium structure of water, we first have to perform a so-called geometry optimization. During the geometry optimization cycle a SCF calculation is performed. This one is used to compute the gradient of the total energy on the nuclei, which is then used

to perform a move of the nuclei. Then the new nuclei positions are feeded into a new SCF calculation. These

steps are repeated until the structure (=positions of the nuclei) does not change anymore.

The input file for a geometry optimization looks like this:

basis=6-31G !basis set

geomtyp=xyz

GEOMETRY=h2o.xyz

hf !Hartree-Fock

optg,savexyz=h2o_opt.xyz

Compare the optimized geometry "h2o_opt.xyz" with your input geometry "h2o.xyz" using VMD. Determine the bond distances and the bending angle in the water molecule. How does it compare to the experimental results? How can we further improve the structure? How many geometry optimization cycles were necessary to converge the geometry? (you could also plot the number of cycles against the corresponding HF energy using "xmgrace" or "gnuplot")

Using the optimized geometry, we are now computing the MOs of water.

Have a look at the input file "h2o_orbitals.com" first and try to analyze what is going on

in this type of calculation. Then type

molpro h2o_orbitals.com

The generated MOs can again be visualized using VMD:

vmd -e orbitals.vmd

Have a look at all occupied MOs and elucidate their role for bonding in the water monomer.

Which atomic orbitals of the individual atoms are predominately participating in each MO?

Compare the MOs with the ones shown on the following web site:

Molecular Orbitals of the Water Monomer

### Caffeine

Finally we are gonna to compute the electronic structure of the well-known

Caffeine molecule, which is a xanthine alkaloid compound that acts as a stimulant in humans.

The structure of this biomolecule is shown to the left.

Type to start the HF caclulation:

molpro caffeine.com &

tail -f caffeine.out

and observe the information that is written to the screen.

How long does one SCF cycle take on average? How many electrons are involved during this computation?

The so-called HOMO (highest occupied molecular orbital) can be visualized using VMD:

vmd -e orbitals.vmd

How many occpuied MO are there?

If there is some time left, try to compute all orbitals by changing "HOMO" to "ALL" in the input file.

Then visualize them along the same lines as in the previous example.