Computational Chemistry/Applications of molecular quantum mechanics

From Wikibooks, open books for an open world
Jump to: navigation, search

Previous chapter - Geometry optimization

Frontier Orbital Theory[edit]

(The Fukui Theory of Reactivity and Selection)

One of the most important pieces of information which can be obtained from the molecular eigenvectors is the Fukui Reactivity Indices. For normal conditions chemical reactivity there are two orbitals which are specially important, the occupied orbital which is highest in energy and the unoccupied orbital which is lowest in energy. These are known as the Frontier Orbitals as they are either side of the occupied / virtual divide. (The empty orbitals which are produced from diagonalising the Fock-Matrix are known as the virtual orbitals.)

The terminology your will come across is the HOMO, Highest Occupied Molecular Orbital, and the LUMO, Lowest Unoccupied Molecular Orbital. Initially nucleophiles attack molecules by placing their surplus electrons, (typically a lone- pair orbital), into the LUMO. The attacking lone pair orbital is then involved as an electron donor and the LUMO an electron acceptor. This is a description in the terminology of Perturbation Molecular Orbital theory. Of course there is a symmetry here in that what is nucleophilic attack by one molecule is electrophilic attach by the other.

Chemical sense is used to decide which label to place on the reaction but we now have a quantitive way to label the reaction in terms of the HOMO and LUMO energies of the reacting species.

In the Fukui theory reactivity indices are calculated from the square of the coefficients of the frontier orbitals. The nucleophilic reactivity indices come from the LUMO coefficients and the electrophilic from the HOMO. For free radical attack the indices are

(c^2_{(HOMO)} + c^2_{(LUMO)}) / 2

Hückel Calculation on Pyridine Output[edit]

  * Pyridine                                                    *

  N C2/C2 C3/C3 C4/C4 C5/C5 C6/C6 N/*
  E(N) = 1.5


      N           C2          C3          C4          C5


  *                            *
  *          6 ATOMS           *
  *                            *
  *          6 ELECTRONS       *
  *                            *


   N    C2   C3   C4   C5   C6

  N        E    1    0    0    0    1

  C2       1    X    1    0    0    0

  C3       0    1    X    1    0    0

  C4       0    0    1    X    1    0

  C5       0    0    0    1    X    1

  C6       1    0    0    0    1    X


         1           2           3           4           5

      2.525982    1.430795    1.000000   -0.594211   -1.000000



  EIGENVECTOR MATRIX                                        PART    1

         1           2           3           4           5

  N          0.754598   -0.432035    0.000000   -0.443012    0.000000
  C2         0.387102    0.014950   -0.500000    0.463881   -0.500000
  C3         0.223215    0.453424   -0.500000    0.167369    0.500000
  C4         0.176735    0.633808    0.000000   -0.563333    0.000000
  C5         0.223215    0.453424    0.500000    0.167369   -0.500000
  C6         0.387102    0.014950    0.500000    0.463881    0.500000

  EIGENVECTOR MATRIX                                        PART    2


  N         -0.218330
  C2         0.367074
  C3        -0.465370
  C4         0.499708
  C5        -0.465370
  C6         0.367074


                  4.46 eV
                277.79 nM
                 430.7 kJ per Mole


      N           C2          C3          C4          C5

     -0.512145    0.199857   -0.010837    0.134105   -0.010837



  N    --- C2           0.5713
  N    --- C6           0.5713
  C2   --- C3           0.6864
  C3   --- C4           0.6537
  C4   --- C5           0.6537
  C5   --- C6           0.6864



      N           C2          C3          C4          C5

     -0.708404    0.234672    0.211150   -0.183239    0.211150



  N    --- C2           0.3658
  N    --- C6           0.3658
  C2   --- C3           0.5140
  C3   --- C4           0.5594
  C4   --- C5           0.5594
  C5   --- C6           0.5140

The reactivity of pyridine is discussed in Streitwieser and Heathcock p1022 onwards. It can be seen from the charge densities that there is positive charge at the ortho (2,5) and para (4) positions, negative charge at the meta (3) position. Electrophiles overwhelmingly attack at the meta position. Both the charge density and the frontier orbitals reinforce here. (The HOMO orbital has equal 0.5 coefficients at ortho and meta positions, but ortho is deactivated by being positively charged.

Nucleophiles attack the LUMO orbital at C2 or C6 (ortho). The coefficients of LUMO at ortho and para are about equal but the ortho position has a greater positive charge which presumably makes it more attractive to electrophiles. (This discussion is about the beginning of the reaction. Sometimes factors which affect the transition state are more important and the predictions from charge densities etc. may not apply when the reaction has proceeded nearer the energy barrier.) In many aromatic systems this kind of frontier orbital analysis works well though. With respect to our discussions of ab-initio quantum chemistry often working in an empty universe with only two molecules involved, it is quite encouraging that Hückel Theory can give rationalizations of reactions which are going on in nitrating mixture, (conc. sulphuric and nitric acids), at 100C plus.

Charge Transfer Interactions[edit]

The so called charge transfer complexes are particularly amenable to frontier orbital analysis. In these complexes the acceptor's LUMO is energetically and geometrically accessible to the donor's HOMO and a small fraction of an electron (typically about 0.05) is passed across. The complex has more conjugation than the separated systems and this often moves the lowest excitation into the visible region so charge-transfer complexes are frequently strongly coloured. ( Streitwieser and Heathcock P.832.

Properties Which can be Calculated by Quantum Mechanics[edit]

Molecular Geometries and Force Constants.[edit]

We have talked about this previously. It is of overriding importance as molecules spend most of their time in the ground state equilibrium geometry, so molecular geometry optimization is usually the first process in any modelling calculation. The force constants lead to the normal coordinates which describe the molecular vibrations. The traversing of a reaction coordinate is sometimes described in rather sophisticated treatments as a molecular vibration of the super-molecule. (When calculating interactions between two molecules the geometry which represents the two molecules in proximity is called the supermolecule. The interaction energy is the energy of the supermolecule minus the sum of the energies of the two free species.)

Atomic Population Analysis and Dipole Moments.[edit]

Atomic population analysis allows the assignment of \delta^+ and \delta^- charges to the component atoms. This is the charge separation caused by the different electronegativities of atoms in different environments. e.g. Azulene which has only carbon and hydrogen atoms in it but has a largeish dipole moment caused by quantum effects. This chemical phenomenon is even accessible to the simplest quantum method, H\"uckel Theory.

As well as the contribution to the dipole from atomic charges there is another contribution caused by distortion of charge on each atom. The semi-empirical programs separate the dipole into these components by labelling them 'POINT-CHARGE' and 'HYBRID' contributions.

Extract from Nitrobenzene MOPAC / AM1 Run[edit]

(Mopac see (Stewart, 1990).)


    1          C          -0.0880          4.0880
    2          C          -0.1397          4.1397
    3          C          -0.0673          4.0673
    4          C          -0.1312          4.1312
    5          C          -0.0678          4.0678
    6          C          -0.1394          4.1394
    7          H           0.1443          0.8557
    8          H           0.1483          0.8517
    9          H           0.1711          0.8289
   10          H           0.1710          0.8290
   11          H           0.1484          0.8516
   12          N           0.5674          4.4326
   13          O          -0.3586          6.3586
   14          O          -0.3585          6.3585
    DIPOLE           X         Y         Z       TOTAL
    POINT-CHG.    -2.598    -4.652    -0.037     5.328
    HYBRID         0.044     0.078     0.001     0.090
    SUM           -2.554    -4.574    -0.036     5.239

Dipoles are often caused by electronegativity differences between the component atoms of the molecule. In some cases however dipoles are caused by entirely electronic quantum effects where the molecule is composed of nothing but carbon and hydrogen. A good example of this is azulene, a fused pair of a five and a seven membered ring. The H\"uckel 4~N~+~2 rule predicts the stability of C6H6 (benzene), C7H7(+) (Tropylium), and C5H5(-) (cyclopentadienyl). The chemistry of these compounds is in no conflict with the theory, C5H5(-) being a stable very common ligand in inorganic chemistry, (it is abbreviated Cp), Ferrocene Fe(C5H5)2 being the first such 'sandwich' compound to be made. Tropylium forms stable ionic salts, e.g. C7H7(+)Br(-).


  C1          C2          C3          C4          C5

  -0.027428   -0.027428    0.145054    0.013553    0.129999

  C6          C7          C8          C9          C10

  0.013553    0.145054   -0.172879   -0.046600   -0.172879

                    .                         .
   ....              3(+)                        ..
      .. 4... ... ...  ...                    ..
        .                 ...                10(-)
       .                     ... 2... ... ... ..
      .                          .              ..
     .                           .                .
  ........5(+)                         .                 ..9........
    ..                           .                 ..
     ..                          .               ..
      ..                         1              .
       ..                    ...  ... ... ... 8(-)
        .6                ...                 .
      ... ... ... ...7(+).                    ..
    ..                .                        ..

Note the alternation of charges with bond connectivity, which
is typical of quantum effects.

  C1   --- C2           0.4009
  C1   --- C7           0.5858
  C1   --- C8           0.5956
  C2   --- C3           0.5858
  C2   --- C10          0.5956
  C3   --- C4           0.6640
  C4   --- C5           0.6389
  C5   --- C6           0.6389
  C6   --- C7           0.6640
  C8   --- C9           0.6560
  C9   --- C10          0.6560

Electronic Spectra[edit]

Electronic spectra can be investigated by several methods the simplest of which of which have the disadvantage that the wavefunction being used is better for the ground state of the molecule than the excited state and this introduces some imbalance into the calculation. Many ab-initio programs now have MCSCF calculations to explore such things analytically. It is however half a PhD project to do it properly.

Ionisation Potentials[edit]

Ionisation potentials are often calculated by the approximate Koopmans' Theorem. This says that the ionisation potential of an electron in a given orbital is just minus the orbital energy. This is simple both to calculate and to interpret but incorrect in that it takes no account of relaxation processes which occur. When ionisation has happened the remaining N-1 electrons can relax to see more of the positive nuclei rather than just staying where they are. This relaxation is effectively instantaneous on the laboratory timescale and so the Koopmans energies are systematically too small. This is the inverse of the process which makes excitation energies in general too large when calculated in the virtual orbital (VO) approximation.

The experimental technique for looking at ionization potentials is photoelectron spectroscopy. This comes in several varieties with different acronyms such as ESCA. The core (1s) ionization potential, despite being tightly bound to the nucleus and taking no part in chemical bonding is still subject to a chemical shift. The simplest example of this is the azide anion, N3. It is charged

   - + -

leading to two photoelectron spectrum peaks at -399eV and -405eV in the ratio 2:1. The practical non-SI unit of spectroscopy at this energy range is the electron volt (eV), 1 eV = 96.4853 kJ per mole, and is about a quarter of the energy needed to break a chemical bond.

Electron Affinities and Excited States[edit]

The same problem of using MO coefficents optimised by the SCF procedure for the ground state configuration applies to the Electron Affinity and to using the virtual orbitals as prototype excited states. Nevertheless this is something which we often do as it is the simplest picture and can be worked on using our readily available ground state semi-empirical wavefunction.

The definitive solution to the above problems is through the use of the electron and polarization propagators.


naphalene and azulene[edit]

Use a molecular gui like Macromodel make trial geometries for naphalene and azulene. Optimise the geometry. Do a single point calculation at the Huckel level and at the AM1 semiempirical level. (A single point calculation is the quantum mechanics at a pre-specified geometry, in this case the MM2* geometry, often either i) the gas phase experimental geometry, ii) the crystal structure geometry from X-ray or neutron diffraction.

Look at the orbital energies and node patterns. Locate both crystal structures on the database and compare their characteristics with the MM2* geometry. Does azulene have C2v symmetry? It should have experimentally in the gas phase but some methodologies break the symmetry and produce only Cs symmetry with alternating single and double bonds. There is a paper on the solution of this problem at the Hartree-Fock level and with correlated methods. (Hartree-Fock gets it wrong.) MP2 and above is right. If you were sufficiently interested you could use a literature search program to find out if the correlated single-determinant method DFT, (Density Functional Theory), gets this structure right.

VSEPR and electronic properties mini project[edit]

This is a mini-project rather than a short problem and involves running the abinitio program of your choice.

The motivation of this computer experiment is to draw together various threads from chemistry about shape, symmetry, electronegativity and electric moments so you will have to consider things like why UF6 is a gas and water and HF hydrogen bond but in the context of the rather dry looking numbers from computer programs which actually mean a lot and can match experiment.

The experiment involves calculation of the electronic structures of CH4, NH3, H2O, HF, BeH2, BF3, H3N-BF3, SF6, XeF4, XeF4O and IF7.

In order to complete this experiment you will need access to a modern molecular quantum mechanics program. You will need to do a bit of literature work and use some common sense to get starting points and experimental numbers for comparison with theory. Starting molecular geometries can always be either calculated or estimated as the ability to guess what the value of a number might be approximately is a neccessary scientific skill. (Crude examples of this are estimating a bond length as the sum of the covalent radii of the atoms, regardless of the environment of the atoms in the current molecule or saying that the chemical shift of an atom is the usual average for that functional group, in the same way as we habitually use average bond energies.) Reference where you get your covalent radii from, there are several sources. You can start with the NH3, BF3 and H3N-BF3 and IF7 (the difficult unusual molecule).

Unless you have a GUI you will need to learn a little bit about a file editor to prepare data files. vi is one of several, pico another. You typically prepare a file for running and submit it to the batch system. (Though many calculations will run in seconds, when preparing work for publication you must always use the best calculation you can afford, and these are typically several hours or even days long.) The quality of the calculation depends on the number of atom based basis orbitals you put in and the cost of the calculations typically increases by the number of basis orbitals to the power 3 and a bit. (You double the basis and the calculation takes 8 times longer.) In a Self Consistent Field (SCF) calculation the usual procedure is to store the many electron repulsion integrals created by the 1/r(1,2) operator in the Hamiltonian on disc. This disc space then becomes a limiting resource which stops the improvement of the calculation. If you want to use lots of f-orbitals on I and d-orbitals on F in IF7 it will fail by filling the disc. Even on NH3 you might use the unfilled d-orbitals in the basis because they allow the nitrogen atom (Another way of looking at this is that the shape of the lone pair and to a lesser extent the tighter bonding orbitals is improved by the extra atomic orbitals.) to be polarized in its bonded environment. There is a mode of execution of most quantum chemistry programs called direct. In this case electron repulsion integrals are calculated on the fly and not stored. This causes the run to take much longer but not to require very much disc space. (This same technique is always used when calculating the 3N-6 position derivatives of electron repulsion integrals during a geometry optimisation as by using a matrix inversion technique these are only required once and not iterated with as in a SCF calculation). If you are really fanatical about doing the definitive calculation on a molecule you might have to use direct and a 32 hour run to use the largest basis to get the best possible calculation. This is not however necessary to get something good out of this experiment.

Estimate the bond lengths in these molecules by adding together the covalent radii. Find the experimental geometries of as many of them as you can. (One of the best compilations of geometries is Landolt-Börnstein. Use your abinitio program to optimize the geometries of the ones you cannot find an experimental geometry for. Here is the experimental geometry of XeF4O below. The units are Angstrom, 10 to the -10 of a metre. You will notice the experimental precision of the Xe=O bond length, and the artificially high precision of the fluorine coordinates which have been subjected to a trigonometrical calculation using 64-bit floating point accuracy! (You will also notice the sloppy way the editor has not recorded the terminal page number of the pre-internet paper and so you might have to go to the library and look it up at some time.)

XeF4=O   J.F. Martins, E.B. Wilson,  J Mol Spectrosc 26, 410 (1968)
#   Xe=O   1.703   Xe-F  1.900  O=Xe-F 91.8degree
O               0.0000000000    0.0000000000    1.703
XE              0.0000000000    0.0000000000    0.0000000000
F               1.8990624647    0.0000000000   -0.0596804422
F              -1.8990624647    0.0000000000   -0.0596804422
F               0.0             1.8990624647   -0.0596804422
F               0.0            -1.8990624647   -0.0596804422

Preparing data for programs like GAUSSIAN, DALTON, GAMESS or CADPAC can be difficult If you have to use the raw data structure. Most programs now have a GUI by which the data preparation can be simplified, even automated.

Answer the following questions in a short report. Give some references and use a style which includes the paper titles and both starting and ending page numbers, (just to get used to the fiddly style used by some journals which waste your time if you loose the end page number!). Many Web and CD journals use this style as they are not wasting paper but it actually does make the reference list easier to navigate.

(Q1) Report the first non-zero moment of these molecules i.e. if the molecule has a dipole moment, the dipole, if it has a no dipole, the quadrupole. Two of these molecules are so symmetrical that they have no dipole or quadrupole. What is the name of the next moment, (rather obscure).

Note that for the dipole moment of XeF4O it is not obvious which direction the dipole will point as both fluorine and oxygen are electronegative and the O-Xe-F angle is close to 90 degrees. However as there is no symmetry along the O=Xe bond so there must be a non-zero dipole.

(Q2) What effect does the lack of a dipole or quadrupole have on the physical properties of these molecules?

(Q3) If there is no dipole moment what forces cause the molecule to condense from gas to a liquid, (or in the case of CO2 straight to a solid at 1 atmosphere). What is the relevance of this to the separation of the isotopes of Uranium? Which two factors fortunately make this extreemly difficult, (hint. somebody's Law is important in one, how you make Uranium a gas in the other)?

The orbital energies are an approximation to the ionization potential. (Technically this is known as Koopmans' Theorem.) The orbital energies are usually in atomic units, but experimentalists use the electron volt (eV).

(Q4) See which electron comes out first in NH3, HCl, IF7 and XeF4. Comment on what you expected and what actually happens.

(Q5) SF6 is 23900 times more potent as a greenhouse gas than CO2. Calculate its infrared spectrum and look at the orbital energies to see how the binding of electrons influences the molecule's reactivity. Examine the vibrational normal modes with respect to the group tables in a book such as Atkins Physical Chemistry or Cotton Group Theory. See how the degeneracies in the group table correspond to both the vibrational degeneracies and the degeneracies in the ionization potentials and therefore the photoelectron spectrum.

Confirm the VSEPR rules by starting BF3 and AlF3 from trigonal geometries like ammonia. Hopefully this optimization will go planar! Calculate the vibrational spectrum of these molecules in their correct geometry.

(Q6) Calculate the vibrational spectrum of flat ammonia and see what happens.

There are methods of estimating the moments of these molecules from electronegativity and shape alone. In general there are two geometries of the molecule we are interested in. The experimental equilibrium geometry and the optimised geometry of the biggest calculation we can afford.

(Q7) Calculate the proton and fluorine shieldings in the molecules which have F and H.

(Q8) Calculate the (H,H) spin-spin coupling constant in NH3 and the (F,F) constant in XeF4, using the SCF wavefunction (as standard). If you have the computing capacity you could calculate the (H,H) spin-spin coupling constant with an expensive electron correlated calculation, which is actually necessary to get good answers for this exotic property.

Some Experimental Geometries (VSEPR)[edit]

NH3 <HNH> 104.9067 theta= 113.7213 R=0.1.024275 G.H.F.D.&A.J.S, J.C.P. 75,1253.
N    0.0            0.000000    0.1449084
H1    0.93773726    0.0           -0.26714530
H2             -0.4688686300    0.8121042892   -0.2671453000
H3             -0.4688686300   -0.8121042892   -0.2671453000

BF3   Kuchitsu K, Konaka S. JCP 45 4342 (1966)
B               0.0000000000    0.0000000000    0.0000000000
F               1.3120000000    0.0000000000    0.0000000000
F              -0.6560000000    1.1362253298    0.0000000000
F              -0.6560000000   -1.1362253298    0.0000000000

AlCl3 2.049 V.P. Spiridonov, A.G. Gershikov, E.Z. Zasorin, N.I. Popenko,
# A.A.Ivanov, L.I. Ermolayeva High. Temp. Sci. 14 (1981) 285.
Al 0 0 0
Cl  2.049
CL             -1.0245000000    1.7744860524    0.0000000000
CL             -1.0245000000   -1.7744860524    0.0000000000

Notice that when dealing with geometries of ammonia and phosphine the most useful angle is the lone pair to X-H angle, (theta), not the H-N-H angle.

Bibliography for VSEPR[edit]

(1) A Chemistry VSEPR website:

(2) Different editions of J. E. Huheey, Inorganic Chemistry.

(3) C. J. Cramer, Essentials of computational chemistry: theories and models, (Chichester, John Wiley, 2002).

(4) F. Jensen, Introduction to Computational Chemistry, (Wiley, Chichester, 1999).

(5) The DALTON website:

(6) Landolt-Börnstein, New series II, Volume 7, Structure Data of Free Polyatomic Molecules, ed. K. H. Hellwege and A. M. Hellwege (Springer, Berlin, 1976). Landolt-Börnstein, New series II, Volume 15 Supplement to Volume II/7, Structure Data of Free Polyatomic Molecules, ed. K. H. Hellwege and A. M. Hellwege (Springer, Berlin, 1987).


  • J. J. P. Stewart, J. Comp. Aided Mol. Design, Vol. 4, 1 (1990).
  • A. Streitwieser and C. H. Heathcock,Introduction to Organic Chemistry, Third Edition, p.236, (Macmillan,New York,1985).

Next Chapter - Applications of molecular modelling