Quantum Chemistry tutorial


Quantum chemistry computations with Molpro (original course material by Dr. U. W. Schmidt)


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
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 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:

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 with a new input file:

basis=6-31G
r=1.0 ang
geometry={angstrom;
          h1;
          h2,h1,r}
hf

{cube,h2
orbital,all
}

and execute molpro again:

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 -f *.cube

will load an existing script that loads the relevant MOs.
To inspect the orbitals, go to Graphics, Representations. Select drawing method CPK. Now you see the two H atoms.
Next, Create Rep (which copies the present representation), and select Drawing method isosurface. With the sliders you can decide the isovalue that is plotted. VMD can show the positive or negative values of the orbital expectation values. To see both, copy the representation Create Rep, and select the negative isovalue points by manualy canging the sign of the Isovalue.

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'
}

h2o

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.

basis=6-31G                  !basis set

geomtyp=xyz
GEOMETRY=h2o_opt.xyz

hf                           !do scf

{cube,h2o
orbital,all
}


Then type

molpro h2o_orbitals.com 

The generated MOs can again be visualized using VMD:

vmd -f *.cube

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

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.




Here are the cartesian coordinates:

  24
HF-SCF000/STO-3G  ENERGY=-667.73565398
C          0.7948047874       -0.2673475369        0.0000000000
C         -0.2830459813       -1.0890658393        0.0000000000
N          1.9143454786       -1.1136247133        0.0000000000
N          0.0798949129       -2.4417816138        0.0000000000
N         -1.6029576700       -0.6005693362        0.0000000000
C          1.4073520392       -2.3930182886        0.0000000000
C          0.7125900943        1.2091306058        0.0000000000
C         -1.8326022679        0.8073510174        0.0000000000
N         -0.6606787360        1.6419581534        0.0000000000
O          1.6414281440        2.0047021485        0.0000000000
O         -2.9562106712        1.2845856383        0.0000000000
C         -0.9238244722        3.0918864981        0.0000000000
C         -2.7666488438       -1.4999441112        0.0000000000
C          3.3316990743       -0.7109956220        0.0000000000
H          2.0612180469       -3.2586428364        0.0000000000
H         -1.4903555584        3.3842737788       -0.8850934873
H         -1.4903555584        3.3842737788        0.8850934873
H          0.0364850129        3.6028877898        0.0000000000
H         -3.3823201570       -1.3336563127       -0.8857759769
H         -3.3823201570       -1.3336563127        0.8857759769
H         -2.4015120487       -2.5249243967        0.0000000000
H          3.9471620995       -1.6091240819        0.0000000000
H          3.5561814227       -0.1173493821        0.8876452774
H          3.5561814227       -0.1173493821       -0.8876452774

and here is what the input file should look like:

 basis=6-31G         ! basis 

 geomtyp=xyz
 GEOMETRY=caffeine.xyz

 set,charge=+0

 gdirect

 hf
 
{cube,caff
orbital,homo
}




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 -f  caff_orbital_10.2.cube

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.

caff