This solver will only work if all \(p\) orbitals are equivalent (it will not consider for instance perturbations due to nitrogen in pyridine). The Huckel theory assumptions are: \[ H_{ij}= \begin{cases} \alpha,\ i=j\\ \beta,\ i \text{ adjacent to } j\\ 0,\ \text{otherwise}\end{cases} \] with scaling such that \(\alpha = 0\) and \(\beta =-1\). As such, it can be seen that the resulting Huckel matrix is the negative of the adjacency matrix \(A\) where: \[ A_{ij}= \begin{cases} 1,\ i \text{ adjacent to } j\\ 0,\ \text{otherwise}\end{cases} \]
generate_smiles (l_or_r:str, n:int)
generate the smiles of either a straight chain or ring polyene, with n atoms. All carbons will be sp2 hybridised. For linear molecules with an odd number of atoms it will return the anion For rings with 4n+1 atoms it will return the anion eg C5H5- For rings with 4n+3 atoms it will return the cation eg C7H7+
| Type | Details | |
|---|---|---|
| l_or_r | str | linear or ring |
| n | int | the number of atoms in the molecule |
Huckel_solve (matrix)
From a Huckel matrix input, solve for the Huckel pi system Returns a dictionary of energy levels with the associated (possibly degenerate) wavefunctions
Huckel (SMILES:str='c1ccccc1')
The solution to the Huckel equation for a molecule given as SMILES
This works in a straightforward manner for chains and rings
mol = generate_smiles('linear', 10)
mol = Huckel(mol)
print(mol)
mol.moleculeHuckel Energies (degeneracy) for C=CC=CC=CC=CC=C: [-1.919 (1)] [-1.683 (1)] [-1.310 (1)] [-0.831 (1)] [-0.285 (1)] [0.285 (1)] [0.831 (1)] [1.310 (1)] [1.683 (1)] [1.919 (1)] mol.plot()
mol = generate_smiles('ring', 6)
mol = Huckel(mol)
print(mol)
mol.moleculeHuckel Energies (degeneracy) for C1=CC=CC=C1: [-2.000 (1)] [-1.000 (2)] [1.000 (2)] [2.000 (1)] mol.plot()
Similarly, it works for more complkex 3d structures (though rdkit struggles to show their structure in the output skeletal diagram)
tetrahedrane_smiles = 'C12C3C1C23'
cubane_smiles = 'C12C3C4C1C5C2C3C45'
dodecahedrane_smiles = 'C12C3C4C5C1C6C7C2C8C3C9C4C1C5C6C2C7C8C9C12'
fullerene_smiles = "c12c3c4c5c2c2c6c7c1c1c8c3c3c9c4c4c%10c5c5c2c2c6c6c%11c7c1c1c7c8c3c3c8c9c4c4c9c%10c5c5c2c2c6c6c%11c1c1c7c3c3c8c4c4c9c5c2c2c6c1c3c42"# rdkit does consider all bonds, even if all are not shown on the output skeletal diagram: this is mostly apparent for cubane
molecule = Chem.MolFromSmiles(cubane_smiles)
mat = -Chem.GetAdjacencyMatrix(molecule)
print(mat)[[ 0 -1 0 -1 0 -1 0 0]
[-1 0 -1 0 0 0 -1 0]
[ 0 -1 0 -1 0 0 0 -1]
[-1 0 -1 0 -1 0 0 0]
[ 0 0 0 -1 0 -1 0 -1]
[-1 0 0 0 -1 0 -1 0]
[ 0 -1 0 0 0 -1 0 -1]
[ 0 0 -1 0 -1 0 -1 0]]tetrahedrane = Huckel(tetrahedrane_smiles)
print(tetrahedrane)
tetrahedrane.plot()Huckel Energies (degeneracy) for C12C3C1C23: [-3.000 (1)] [1.000 (3)] 
cubane = Huckel(cubane_smiles)
print(cubane)
cubane.plot()Huckel Energies (degeneracy) for C12C3C4C1C5C2C3C45: [-3.000 (1)] [-1.000 (3)] [1.000 (3)] [3.000 (1)] 
dodecahedrane = Huckel(dodecahedrane_smiles)
print(dodecahedrane)
dodecahedrane.plot()Huckel Energies (degeneracy) for C12C3C4C5C1C6C7C2C8C3C9C4C1C5C6C2C7C8C9C12: [-3.000 (1)] [-2.236 (3)] [-1.000 (5)] [0.000 (4)] [2.000 (4)] [2.236 (3)] 
fullerene = Huckel(fullerene_smiles)
print(fullerene)
fullerene.plot()Huckel Energies (degeneracy) for c12c3c4c5c2c2c6c7c1c1c8c3c3c9c4c4c%10c5c5c2c2c6c6c%11c7c1c1c7c8c3c3c8c9c4c4c9c%10c5c5c2c2c6c6c%11c1c1c7c3c3c8c4c4c9c5c2c2c6c1c3c42: [-3.000 (1)] [-2.757 (3)] [-2.303 (5)] [-1.820 (3)] [-1.562 (4)] [-1.000 (9)] [-0.618 (5)] [0.139 (3)] [0.382 (3)] [1.303 (5)] [1.438 (3)] [1.618 (5)] [2.000 (4)] [2.562 (4)] [2.618 (3)] 
It is also possible to manually overwrite the Huckel matrix in order to deal with species such as pyridine:
pyridine = Huckel('n1ccccc1')
pyridine.moleculeThe Huckel matrix at the moment is still only the negative of the adjacency matrix, and so is the same as that of benzene
pyridine.matrixarray([[ 0, -1, 0, 0, 0, -1],
[-1, 0, -1, 0, 0, 0],
[ 0, -1, 0, -1, 0, 0],
[ 0, 0, -1, 0, -1, 0],
[ 0, 0, 0, -1, 0, -1],
[-1, 0, 0, 0, -1, 0]], dtype=int32)We can overwrite this matrix: taking the assumption in the A4 course, that \(\alpha_{N} = \alpha+\beta/2\) and \(\beta_{NC}=0.8\beta\)
pyridine.matrix = np.array([[ -1/2, -0.8, 0, 0, 0, -0.8],
[-0.8, 0, -1, 0, 0, 0],
[ 0, -1, 0, -1, 0, 0],
[ 0, 0, -1, 0, -1, 0],
[ 0, 0, 0, -1, 0, -1],
[-0.8, 0, 0, 0, -1, 0]])print(pyridine)
pyridine.plot()Huckel Energies (degeneracy) for n1ccccc1: [-1.954 (1)] [-1.062 (1)] [-1.000 (1)] [0.667 (1)] [1.000 (1)] [1.849 (1)] 