Previous

3. Natural Molecular Orbital Analysis

The occupied canonical orbitals that are obtained from a Hartree-Fock (HF) calculation are generally delocalized. In many situations, this is undesirable because it is difficult to attach a chemical interpretation to these molecular orbitals. However, the set of canonical molecular orbitals can be transformed into equally valid set of localized HF molecular orbitals by a (unitary) transformation that preserves orthonormality. These localized molecular orbitals (LMOs) correspond to chemically familiar concepts including the core orbitals on the heavy atoms, bond orbitals, and the lone pair orbitals. The latter can be either sp hybrids or σ and π orbitals. The pair of images below illustrates the canonical (left) and natural (left) highest occupied molecular orbital (HOMO) for acetone.

Acetone: Canonical HOMO Acetone: Natural HOMO

The localization is usually carried out as part of a so-called Natural Bond Orbital (NBO) analysis. Natural Bond Orbitals can be obtained by diagonalization of localized blocks of the density matrix. They correspond to the Lewis structure representation of the molecule describing "localized electrons" in bonds and lone pairs. For example, notice that the HOMO of acetone can be clearly interpreted as the orbital holding one of the lone pairs.

If resonance is important the molecule cannot be represented by a single Lewis structure and the NBO analysis has to be extended to include combinations of the relevant NBOs to form Natural Localized Molecular Orbitals (NLMOs). The Gaussian computer code automatically makes this extension. However, for molecules where there is no significant resonance, the NLMOs and NBOs are essentially the same. For example, the NLMOs of methane consist of one carbon core orbital and four equivalent sp3 hybrid C-H bond orbitals.

Gaussian performs the NBO/NLMO analysis via the keyword Pop=NBO. In order to visualize the NLMOs with MOLDEN, a two-step job must be constructed. The first step performs the NBO/NLMO analysis and saves the NLMOs into the checkpoint file using the keywords Pop=(NBORead,SaveNLMOs). The NBORead directive requires an extra line after the Z-matrix that in this case reads $NBO AONLMO $END. For example, the input file to generate NLMOs for methane might look like this:

%Mem=32MW
%Chk=ch4.chk
#P HF/3-21G Opt GFINPUT IOP(6/7=3) Pop=(NBORead,SaveNLMOs)

Methane: NMO Analysis; write NLMOs to Checkpoint file

0  1
 C
 H     1     hc
 H     1     hc        2     hch
 H     1     hc        2     hch       3     dih       0
 H     1     hc        2     hch       3     -dih      0

  hc                    1.08618
  hch                 109.47122
  dih                 120.

$NBO AONLMO $END

The second step reads the NLMOs from the checkpoint file and writes them into a Gaussian output file in a format suitable for visualization with MOLDEN:

%Chk=ch4.chk
#P HF/3-21G Geom=Check IOP(6/7=3) GFINPUT Guess=(Read,Only)

Methane NLMO: Get NLMOs from Chk to Out file

0  1

Further information on NBO:

4. Wavefunction instability

Near the equilibrium geometry most molecules that we encounter in everyday life are well-approximated by a restricted Hartree-Fock (RHF) wavefunction. However, for some molecules near the equilibrium, such as molecular oxygen (a diradical), an unrestricted Hartree-Fock (UHF) wavefunction will yield a better approximation. The UHF method is often important for the description of bond-breaking. It is often difficult to know a priori when this circumstance will occur. In order to do so one can carry out an RHF calculation and, then, test to see if the energy can be lowered by removing the double-occupancy restriction, i.e. by using a UHF wavefunction. If that is the case, we say that the RHF wavefunction found by is unstable. The stability can be tested in Gaussian using the Stable keyword. For example, performing a stability analysis for a system consisting of 2 hydrogen atoms 10 Angstroms apart will yield:

 
 ***********************************************************************
 Stability analysis using  singles matrix:
 ***********************************************************************

 Eigenvectors of the stability matrix:

 Eigenvector   1:   Triplet-SGU  Eigenvalue=-0.7216882
       1 ->  2         0.70711

 Eigenvector   2:   Singlet-SGU  Eigenvalue= 0.7216882
       1 ->  2         0.70711
 The wavefunction has an RHF -> UHF instability.  
 
      Molecular Orbital Coefficients
                           1         2
                       (SGG)--O  (SGU)--V
     EIGENVALUES --    -0.10574  -0.05282
   1 1   H  1S          0.70711   0.70711
   2 2   H  1S          0.70711  -0.70711

Notice that the occupied MO coefficients correspond to an orbital that is delocalized over both H atoms whereas we know that, at long internuclear distances there should be one electron localized on each H atom. If the RHF wavefunction is unstable, a lower energy solution may be found by invoking the calculation with the Stable=Opt keyword. If instability is found the RHF wavefunction is allowed to become UHF. Molecular hydrogen is a simple case because it has the same symmetry at all internuclear distances. For an arbitrary molecule the program will also check for instability with respect to reducing molecular symmetry.


 ***********************************************************************
 Stability analysis using  singles matrix:
 ***********************************************************************

 Eigenvectors of the stability matrix:

 Eigenvector   1:   Triplet-SGU  Eigenvalue=-0.7216882
       1 ->  2         0.70711

 Eigenvector   2:   Singlet-SGU  Eigenvalue= 0.7216882
       1 ->  2         0.70711
 The wavefunction has an RHF -> UHF instability.
 Requested convergence on RMS density matrix=1.00D-08 within 128 cycles.
 Requested convergence on MAX density matrix=1.00D-06.
 Requested convergence on             energy=1.00D-06.
 Reduce DIIS space if new energy is above lowest so far.
 Dynamic level shift is off after energy rises.
 Keep R1 and R2 integrals in memory in canonical form, NReq=      818934.
 SCF Done:  E(UHF) = -0.933163700769     A.U. after    7 cycles
	    
      Alpha Molecular Orbital Coefficients
                           1         2
                           O         V
     EIGENVALUES --    -0.46658   0.30802
   1 1   H  1S          0.00000   1.00000
   2 2   H  1S          1.00000   0.00000

Notice that now the MO coefficients show a physically correct behavior. You can retrieve the correct MO coefficient into a new file using the Guess=(Read,Only) keyword.


Previous

Materials by Dr. Kalju Kahn and Bernie Kirtman, Department of Chemistry and Biochemistry, UC Santa Barbara. 2005-2012