Roberto Di Remigio, Luca Frediani
With the Hückel method we aim at calculating the energy levels of conjugated $\pi$ systems. The method is limited to planar, conjugated hydrocarbons and includes only the $\pi$ electrons in the system in the quantum mechanical treatment. The $\sigma$ electrons are ignored, as they are at much lower energies and will not affect molecular properties significantly. The molecular orbitals (MOs) are built as a linear combination of atomic orbitals (LCAO), where the atomic orbitals (AOs) in question are the $p_z$ orbitals on each of the carbon atoms.
If $N$ is the number of carbon atoms in the conjugated hydrocarbon, the trial wavefunction is thus written as: \begin{equation} \Psi = \sum_{i=1}^{N}c_i\chi_i \end{equation} where the index $i$ runs over carbon atoms and $\chi_i$ is the $p_z$ orbital on the $i$-th atom in the molecule. The energy, i.e. the expectation value of the many-electron Hamiltonian, is: \begin{equation} E(\mathbf{c}^\dagger, \mathbf{c}) = \frac{\langle \Psi | H | \Psi \rangle}{\langle \Psi | \Psi \rangle} = \frac{\sum_{i,j=1}^{N} c_i^*H_{ij}c_j }{\sum_{i, j=1}^{N}c_i^*S_{ij}c_j} = \frac{\mathbf{c}^\dagger\mathbf{H}\mathbf{c}}{\mathbf{c}^\dagger\mathbf{S}\mathbf{c}} \end{equation} where the Hamiltonian and overlap matrices were introduced: \begin{alignat}{2} H_{ij} = \langle \chi_i | H |\chi_j\rangle \quad& S_{ij} = \langle \chi_i | \chi_j \rangle \end{alignat} The variation principle guarantees that minimization of the energy function will lead to the optimal value of the coefficients, in the chosen LCAO basis. This reduces to the generalized eigenvalue problem: \begin{equation} \mathbf{H}\mathbf{c} = \mathbf{S}\mathbf{c}E \end{equation} where the eigenvector $\mathbf{c}$ is the set of coefficients in the linear combination describing the MO with energy eigenvalue $E$. Finding the eigenvalues can be achieved by calculating the determinant of $\mathbf{H} - E\mathbf{S}$ and setting it to zero: \begin{equation} \det(\mathbf{H} - E\mathbf{S}) = 0 \end{equation}
As already mentions, the molecular orbitals (MOs) are built as a linear combination of atomic orbitals (LCAO), where the atomic orbitals (AOs) in question are the $p_z$ orbitals on each of the carbon atoms. A schematic view of the model is given in the Figure below.
The Hamiltonian and overlap matrices are constructed from a simple set of rules:
Once the connectivity and the values of the $\alpha$ and $\beta$ parameters are known one can set up the Hückel matrix and the correspoding secular problem.
Ethene has two carbon atoms. The Hückel matrix is: \begin{equation} \begin{pmatrix} \alpha & \beta \\ \beta & \alpha \end{pmatrix} \end{equation} and its secular problem: \begin{equation} \det \begin{pmatrix} (\alpha - E) & \beta \\ \beta & (\alpha - E) \end{pmatrix} = (\alpha - E)^2 - \beta^2 = 0 \end{equation} The solution of the second order equation gives the energies for the highest occupied molecular orbital $E_-$ (HOMO) and the lowest unoccupied molecular orbital $E_+$ (LUMO): \begin{equation} E_- = \alpha -\beta \quad E_+ = \alpha + \beta \end{equation} Inserting these energies in the eigenvalue equation will lead to two different systems of equations whose solution determines the coefficient in the MOs linear combinations. Notice that the coefficients are determined up to normalization, meaning that only by imposing a normalization condition the solution to the eigenvalue problem will be unique.
Write a Python function to calculate eigenvalues and eigenvectors for the Hückel matrix of a generic conjugated hydrocarbon. The function should take four parameters:
The connectivity array describes which carbon atoms are connected via a double bond. This is essential, since it determines which entries in the Hückel matrix are nonzero!
As an example, for the ethene molecule $[1, 2]$ would be the connectivity array. For the butadiene molecule one would instead have $[1, 2, 2, 3, 3, 4]$ which means that atom $1$ is connected to atom $2$ by a double bond, atom $2$ and $3$ are connected by a single bond and atom $3$ and $4$ are again connected by a double bond.
Why so? Because butadiene has 4 carbon atoms, hence only two double bonds can be present. Since these have to be alternating carbon $2$ and carbon $3$ must be connected by a single bond.
To help with determining the connectivity array, draw the molecule and number all carbon atoms.
The function will, of course, use NumPy for manipulating matrices and vectors.
The module linalg
includes the functions necessary for linear algebra manipulations. In our case, we're interested in finding eigenvalues and eigenvectors of the Hückel matrix with the eig
function
def huckel(n, alpha, beta, connectivity):
""" Huckel method
n -- number of carbon atoms
alpha -- Coulomb integral
beta -- resonance integral
connectivity -- an array of integers defining the connectivity
"""
import numpy as np
from numpy import linalg as LA
H = np.zeros((n, n))
# Fill the matrix!
## Fill diagonal
## Fill off-diagonal
# Print H
print('Huckel matrix is {}'.format(H))
# Diagonalize H
energies, orbitals = LA.eig(H)
# Sort eigenvalues and eigenvectors in ascending order
sort_perm = energies.argsort()
energies.sort() # <-- This sorts the list in place.
orbitals = evecs[:, sort_perm]
# Return energies and orbitals
return energies, orbitals
In [10]:
def huckel(n, alpha, beta, connectivity):
""" Huckel method
n -- number of carbon atoms
alpha -- Coulomb integral
beta -- resonance integral
connectivity -- an array of integers defining the connectivity
"""
import numpy as np
from numpy import linalg as LA
H = np.zeros((n, n))
# Fill the matrix!
## Fill diagonal
## Fill off-diagonal
# Print H
print('Huckel matrix is \n{}'.format(H))
# Diagonalize H
energies, evecs = LA.eig(H)
# Sort eigenvalues and eigenvectors in ascending order:q
sort_perm = energies.argsort()
energies.sort() # <-- This sorts the list in place.
orbitals = evecs[:, sort_perm]
# Return energies and orbitals
return energies, orbitals
energies, orbitals = huckel(4, -5.0, -75.0, [1, 2, 2, 3, 3, 4])
print('Huckel model energies \n{}'.format(energies))
print('Huckel model orbitals \n{}'.format(orbitals))
We will use the packages chemview
and chemlab
to visualize the molecular orbitals, as in this example. The coefficients are calculated from our small function for the Hückel method. The technicalities of the definition of the basis functions are not important in this context. You can find more details here.
In [1]:
import numpy as np
from chemlab.db.cirdb import CirDB
from chemlab.qc import molecular_orbital
from chemview import enable_notebook, MolecularViewer
enable_notebook()
# To get molecule by name
ethene = CirDB().get("molecule", "ethene")
# Huckel's method results for ethene
orbitals = np.array([[ 0.70710678, 0.70710678],
[-0.70710678, 0.70710678]])
# Initialize molecular viewer
mv = MolecularViewer(ethene.r_array, { 'atom_types': ethene.type_array,
'bonds': ethene.bonds })
# Define rendering type
mv.wireframe()
# Define basis set
# p-functions for carbon from STO-3G basis set
sto3g = [
# atom 1
[
('P', [
(2.941249, 0.155916), (0.683483, 0.607684), (0.222290, 0.391957), ]),
],
# atom 2
[
('P', [
(2.941249, 0.155916), (0.683483, 0.607684), (0.222290, 0.391957), ]),
]
]
# Define molecular orbital to be plotted
# Notice that we pass only the coordinates of the carbon atoms!
coords = np.array([ethene.r_array[0], ethene.r_array[1]])
# We need to fudge the coefficient arrays with zeros (p_x and p_y are absent!)
HOMO_coeff = np.array([0.0, 0.0, orbitals[0][0], 0.0, 0.0, orbitals[0][1]])
LUMO_coeff = np.array([0.0, 0.0, orbitals[1][0], 0.0, 0.0, orbitals[1][1]])
HOMO = molecular_orbital(coords, HOMO_coeff, sto3g)
LUMO = molecular_orbital(coords, LUMO_coeff, sto3g)
# Add isosurface of MO with positive isovalue
mv.add_isosurface(HOMO, isolevel=0.3, color=0xff0000)
# Add isosurface of MO with negative isovalue
mv.add_isosurface(HOMO, isolevel=-0.3, color=0x0000ff)
# This is necessary to actually get a picture of the molecule
mv
In [ ]: