Past semester we revisited the Hückel molecular orbital theory at class. One day I was trying to solve some problems with SymPy, a Python module that is a computer algebra system (CAS), and noticed something curious while getting symbolic solutions for the pentafulvene molecule, as shown below:
Python 3.8.5 (default, Sep 4 2020, 07:30:14)
Type "copyright", "credits" or "license" for more information.
IPython 7.19.0 -- An enhanced Interactive Python.
from sympy import Symbol, Matrix, roots
x = Symbol('x')
M_pentafulvene = Matrix([[x, 1, 0, 0, 0, 0], [1, x, 1, 0, 0, 1], [0, 1, x, 1, 0, 0], [0, 0, 1, x, 1, 0], [0, 0, 0, 1, x, 1], [0, 1, 0, 0, 1, x]])
M_pentafulvene
Out[4]:
Matrix([
[x, 1, 0, 0, 0, 0],
[1, x, 1, 0, 0, 1],
[0, 1, x, 1, 0, 0],
[0, 0, 1, x, 1, 0],
[0, 0, 0, 1, x, 1],
[0, 1, 0, 0, 1, x]])
polynomial = M_pentafulvene.det()
polynomial
Out[6]:
x**6 - 6*x**4 + 8*x**2 + 2*x - 1
roots(polynomial)
Out[7]:
{-1: 1,
1/2 - sqrt(5)/2: 1,
1/2 + sqrt(5)/2: 1,
-(27/2 + 3*sqrt(687)*I/2)**(1/3)/3 - 4/(27/2 + 3*sqrt(687)*I/2)**(1/3): 1,
-(-1/2 + sqrt(3)*I/2)*(27/2 + 3*sqrt(687)*I/2)**(1/3)/3 - 4/((-1/2 + sqrt(3)*I/2)*(27/2 + 3*sqrt(687)*I/2)**(1/3)): 1,
-4/((-1/2 - sqrt(3)*I/2)*(27/2 + 3*sqrt(687)*I/2)**(1/3)) - (-1/2 - sqrt(3)*I/2)*(27/2 + 3*sqrt(687)*I/2)**(1/3)/3: 1}
As you can see, the expressions for three out of five solutions have $I = \sqrt -1$, the imaginary unit, as part of them. Seeing that, initially I was puzzled, as imaginary solutions would make no sense, and thought I made a mistake somewhere, but soon I perceived the imaginary part pratically vanishes when we evaluate the expressions:
I = (-1)**0.5
sqrt = lambda x: x**0.5
for solution in roots(polynomial):
print(eval(str(solution)))
-1
-0.6180339887498949
1.618033988749895
(-2.114907541476756+1.6653345369377348e-16j)
(1.8608058531117035+3.3306690738754696e-16j)
(0.2541016883650525-4.440892098500626e-16j)
There is a tiny imaginary remnant, on the order of $10^{-16}$. I don't know the implementation details of complex numbers in Python, but I suspect this tiny remnant is due some sort of rounding error.
Anyway, as is not that hard to came with polynomials containing integer coeficients that yield imaginary solutions, the trivial example being $x^2 + 1 = 0$. That made me wonder if one could "break" Hückel method by finding a example of a plausible compound whose associated secular determinant results in complex roots with a non negligible imaginary part. I actually tried to brute-force an answer, by writing a script to try solve every determinant for "valid" matrices up to a 7 carbon scaffold. By valid I mean, matrices that are symmetric around the main diagonal, whose lines and columns sum to a number greater than zero and equal or smaller than 3 (to avoid carbons with too many neighbors or none), and that represented fully connected structures:
import copy
import time
from sympy import Matrix, eye, symbols, nroots
def load_symbol(L):
"""
Takes all-numeric square matrix and replaces every entry in the main diagonal
with variable x (SymPy symbol), returning a symbolic matrix.
"""
x = symbols('x')
for i in range(len(L)):
L[i][i] = x
return L
def valid_columns_rows(L):
"""
Takes square matrix (list of n integer lists of size n), and verify if it
could represent a valid Hückel secular determinant, by checking, if the sum
of each row is > 0 and <= 3, then doing the same check to each column, as
each sp2 carbon should be bonded to at least another carbon, and at most to
other 3 coplanar carbons. Returns True if both checks are True, otherwise
returns false.
"""
row_sums = [sum(row) for row in L]
if min(row_sums) > 0 and max(row_sums) <= 3:
column_sums = []
for i in range(len(L)):
ith_column = [row[i] for row in L]
column_sums.append(sum(ith_column))
return min(column_sums) > 0 and max(column_sums) <= 3
else:
return False
def connected(L):
"""
Checks if matrix represents fully connected structure, instead of a set
of disjoined smaller fragments. Returns True if connected, otherwise
returns False. As the Huckël secular determinant matrix has a similar form
to the adjacency matrix of a graph, I used the power sum of adjacency
matrix method given in Tan T Z, Gao S X, Yang W G. Determining the
connectedness of an undirected graph. Journal of University of Chinese
Academy of Sciences, 2018, 35(5): 582-588, available at
http://html.rhhz.net/ZGKXYDXXB/20180502.htm
It consists in calculating:
S = I + M + M^2 + ... + M^(n-1),
where M is the adjacency matrix and I the identity matrix of same size.
If the first row of S has all elements > 0, then the system is connected.
"""
n = len(L)
M = Matrix(L)
S = eye(n)
for i in range(n-1):
S += M**(i+1)
return min(S.row(0)) > 0
def diagonal_slice_generator(n):
"""
Generates list of possible diagonal slices of size n, for a huckel secular
determinant matrix. As after substituting (alpha - E)/beta = x, all elements
off the principal diagonal become either 0 or 1, it just generates every
binary number from 0 to 2**n - 1, pads with zeroes, converts to string
then to integer lists.
"""
string_list = [list(str(bin(i))[2:].zfill(n)) for i in range(2**n)]
integer_nested_list = [[int(char) for char in string] for string in string_list]
return integer_nested_list
def secular2roots(M):
"""
Given a symbolic matrix, tries to solve the Hückel secular determinant,
numerically. If successful, returns solutions as list. If it fails to
converge, returns empty list, and prints a warning.
"""
symbolic_matrix = Matrix(M)
polynomial = symbolic_matrix.det()
try:
D = nroots(polynomial, maxsteps=500)
return D
except:
print('Roots not found for polynomial', polynomial, 'from matrix:', M)
return []
def symmetric_matrix_generator(n, L):
"""
Takes integer n and list L, that must be composed of n sublists containing
only 0 as element. Then grows recursively this list into a square matrix of
size n, by adding pairs of identical diagonal slices above and below the
main diagonal. The main diagonal is composed of the zeroes in the original
list. It's a generator, yielding every possible n x n matrix symmetric
around the main diagonal, until the search space is exhausted. The matrix
must be symmetric because if carbon i is bound to carbon j, then carbon j
must be bound to carbon i too, as chemical bonds are bidirectional.
"""
diagonals = diagonal_slice_generator(n-1)
backup = copy.deepcopy(L)
for cut in diagonals:
for j in range(len(cut)):
L[j+(len(L) - len(cut))].insert(0, cut[j])
L[j].append(cut[j])
if n > 2:
yield from symmetric_matrix_generator(n-1, L)
else:
yield L
L = copy.deepcopy(backup)
def valid_matrix_generator(n):
"""
Takes n as argument, then calls the previous symmetric_matrix_generator
to get candidate square matrices of size n x n. Reject the ones unlikely
to represent valid Hückel secular determinants, due to either not being
fully connected, or to carbon atoms with none or more than 3 neighbors,
yielding the rest.
"""
L = [[0] for i in range(n)]
for e in symmetric_matrix_generator(n, L):
if valid_columns_rows(e):
if connected(e):
yield e
n = 7 #Upper bound on size of matrix (lower bound is 2)
maximum = 0 #Stores largest imaginary part found
counter = 0 #Counts number of valid matrices processed
start = time.time()
for i in range(2, n+1):
for M in valid_matrix_generator(i):
s = load_symbol(M)
numeric_roots = secular2roots(s)
if numeric_roots == []:
pass
else:
if len(numeric_roots) < i:
print('Missing roots in', s, 'found only', numeric_roots)
candidate = abs(max([complex(n).imag for n in numeric_roots]))
if candidate > maximum:
maximum = candidate
print('So far, largest imaginary part found in', numeric_roots, 'from', s)
counter += 1
if counter % 10000 == 0:
print('So far', counter, 'valid matrices processed in', round(time.time() - start), 'sec.')
print('Search space exhausted for n=', i, 'at', round(time.time() - start, 2), 'sec.')
print('Total of', counter, 'matrices processed in', round(time.time() - start), 2, 'sec.')
print('Maximum absolute imaginary part found in any solution set:\n', maximum)
Running this, I got the following output after nearly five hours:
Search space exhausted for n= 2 at 0.0 sec.
Search space exhausted for n= 3 at 0.02 sec.
Roots not found for polynomial x**4 - 6*x**2 + 8*x - 3 from matrix: [[x, 1, 1, 1], [1, x, 1, 1], [1, 1, x, 1], [1, 1, 1, x]]
Search space exhausted for n= 4 at 0.56 sec.
Search space exhausted for n= 5 at 10.17 sec.
So far, largest imaginary part found in [-3.00000000000000, -3.04533612508806e-16, 2.54155613029339e-16, 3.00000000000000, -2.78240672724667e-16*I, 3.33228889079154e-16*I] from [[x, 0, 0, 1, 1, 1], [0, x, 0, 1, 1, 1], [0, 0, x, 1, 1, 1], [1, 1, 1, x, 0, 0], [1, 1, 1, 0, x, 0], [1, 1, 1, 0, 0, x]]
Search space exhausted for n= 6 at 405.84 sec.
So far, largest imaginary part found in [-2.76155718183189, -1.36332823793268, -1.00000000000000, 1.00000000000000, 2.12488541976457, 0.999999999999903 + 1.674006560727e-13*I, 1.0000000000001 - 1.69100888096705e-13*I] from [[x, 0, 0, 1, 0, 1, 1], [0, x, 0, 0, 1, 1, 0], [0, 0, x, 0, 1, 0, 1], [1, 0, 0, x, 0, 1, 1], [0, 1, 1, 0, x, 0, 0], [1, 1, 0, 1, 0, x, 0], [1, 0, 1, 1, 0, 0, x]]
So far 10000 valid matrices processed in 512 sec.
So far 20000 valid matrices processed in 1344 sec.
So far 30000 valid matrices processed in 2142 sec.
So far 40000 valid matrices processed in 3014 sec.
So far 50000 valid matrices processed in 3820 sec.
So far 60000 valid matrices processed in 4678 sec.
So far 70000 valid matrices processed in 5646 sec.
So far 80000 valid matrices processed in 6648 sec.
So far 90000 valid matrices processed in 7686 sec.
So far 100000 valid matrices processed in 8770 sec.
So far 110000 valid matrices processed in 9928 sec.
So far 120000 valid matrices processed in 11113 sec.
So far 130000 valid matrices processed in 12374 sec.
So far 140000 valid matrices processed in 13698 sec.
So far 150000 valid matrices processed in 14996 sec.
Search space exhausted for n= 7 at 16231.49 sec.
Total of 158985 matrices processed in 16231 2 sec.
Maximum absolute imaginary part found in any solution set:
1.674006560727002e-13
After processing more than 150k matrices (the number of actual compounds is smaller due to some of the them corresponding to different labelling on the carbons), the numeric root finder of SymPy failed to find a solution in a single case. I tried to draw it, and I don't believe it represents a valid compound to apply Hückel method, as it probably can't be planar. The closest thing would be a 3D carbon tetrahedron where each face looks like a cyclopropene:
It seems my valid_matrix_generator() doesn't generate only valid examples. I think the problem is the lack of criteria to filter out compounds that can exist only in 3D. It's curious that the numeric solver failed to converge for this fourth degree equation, even if it could find solutions for all the others up to seventh degree. The symbolic solver show that it actually has solutions 1 and -3, corresponding to $(x-1)^3(x+3) = 0$.
Anyway, it failed in finding a single example with solutions containing a non negligible imaginary part. The largest found was of the order of $10^{-13}$. While the value is about 1000x larger than what I found in the pentafulvene example, that's still very small. I think it's another compound that would not actually be planar, and the equivalent 3D structure would correspond to something like this:
In the end, as expected, it looks like it's not possible to get imaginary solutions for a Hückel secular determinant, but I'm not sure. If someone could explain me what mathematically precludes a solution being complex, I would be grateful.
