6
$\begingroup$

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:

enter image description here

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:

enter image description here

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.

$\endgroup$
1
  • 1
    $\begingroup$ Just wanted to say that those are excellent diagrams! $\endgroup$ – taciteloquence Jan 19 at 16:26
7
$\begingroup$

In the Huckel method, you are just generating a very simplified version of the molecular Hamiltonian and determining it's eigenvalues. The molecular Hamiltonian will always be a Hermitian matrix (for the Huckel method, we can be more specific and say it's a real symmetric matrix by construction). Hermitian matrices are guaranteed to have all real eigenvalues (a simple proof is available on Math SE and tends to be covered in early quantum mechanics courses), so you won't be able to find a molecule where the Huckel method would produce imaginary/complex energies.

As for the small imaginary values you are seeing, these are a mainly a consequence of the finite precision offered by floating point numbers. ~$10^{-16}$ is roughly the machine epsilon for double precision arithmetic, so it's not surprising to see some of your results have this small an imaginary part. As for the max value you saw of ~$10^{-13}$, the slightly larger imprecision can be attributed to the various solver routines you use having their own inherent imprecision when working with floating point numbers, leaving you with numerical noise in the imaginary part.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.