A simple tool to create pi systems for both linear and cyclic conjugated alkenes.

hueckel_orbitals.py

#
# Version 0.3.
#   - using scipy.linalg.eigh instead of numpy.linalg.eigh provides more accurate eigenvectors
#   - format now supports "pretty" for linear molecules
#
# Call from command line
#   hueckel_orbitals 6 (creates the pi system belonging to hexatriene)
#   hueckel_orbitals 6 --cyclic (creates the pi system for benzene)
# Further options:
#   --reverse       Reverses the order of the output (lowest energy first instead of highest first)
#   --format=csv    Prints the orbitals line-per-line as a csv file containing relative energy and orbital coefficients

from __future__ import (absolute_import, division, print_function,
                        unicode_literals)
from collections import OrderedDict

import argparse
import numpy as np
from scipy import linalg as LA

ROUND = 5
ETA = 10**(-ROUND)

p_0 = "."
p_P = "o"
p_N = "x"

class piSystem():
    _n = 0
    _cyclic = False
    _data = None

    def __init__(self, N, cyclic = False):
        self._n = N
        self._cyclic = cyclic

    def calculateOrbitals(self):
        # Create Matrix
        M = self.createLinearMatrix() if self._cyclic == False else \
            self.createCyclicMatrix()
        # Get eigenvalues and eigenvectors
        eVal, eVec = LA.eigh(M)

        data = {}
        for i in range(0, self._n):
            val = np.round(-eVal[i], ROUND)
            if val not in data:
                data[val] = []
            data[val].append(eVec[:,i])
        data = OrderedDict(sorted(data.items(), key=lambda t: t[0]))

        self._data = data

    def createLinearMatrix(self):
        M = np.zeros((self._n,)*2, dtype=np.float64)
        for i in range(0, self._n):
            if i-1 >= 0:
                M[i,i-1] = 1
            if i+1 < args.n:
                M[i,i+1] = 1
        return M

    def createCyclicMatrix(self):
        M = np.zeros((self._n,)*2, dtype=np.float64)
        for i in range(0, self._n):
            M[i,i-1] = 1
            if i+1 < self._n:
                M[i,i+1] = 1
            else:
                M[i,i+1-self._n] = 1
        return M

    def printSimple(self, reverse=True):
        self.calculateOrbitals()
        lines = []
        for val in self._data:
            line = []
            for o in self._data[val]:
                o*=-1 if o[-1] < 0 else 1
                for pos in range(0, len(o)):
                    if abs(o[-pos-1]) < ETA:
                        line.append(p_0)
                    elif o[-pos-1] > 0:
                        line.append(p_P)
                    else:
                        line.append(p_N)
                line.append(", ")

            lines.append(''.join(line[:-1]))
        for i in range(len(lines)):
            if reverse:
                print(lines[i])
            else:
                print(lines[-i-1])

    def printPretty(self, reverse=True):
        self.calculateOrbitals()
        items = list(self._data.items())
        items.reverse()
        self._data = OrderedDict(items)

        if self._cyclic == False:
            resolution = 3
            empty = [" "*10 + "|"]
            lines = [
                [" "*8 + "E" + " ^"],
            ]

            # Get smallest delta
            smallest_delta = None
            value_before = None
            for val in self._data:
                if value_before is None:
                    value_before = val
                    continue
                delta = abs(value_before - val)
                if smallest_delta is None or delta < smallest_delta:
                    smallest_delta = delta

            # Prepate output
            value_before = None
            for val in self._data:
                valform = "{0:>9.4} + ".format(val)
                # Create orbitals
                line = []
                o = self._data[val][0]
                for pos in range(0, len(o)):
                    if abs(o[-pos-1]) < ETA:
                        line.append(p_0)
                    elif o[-pos-1] > 0:
                        line.append(p_P)
                    else:
                        line.append(p_N)
                line = "".join(line)
                # Create line (and correct line position)
                if value_before is None:
                    lines.append(valform + line)
                else:
                    delta = round(abs(value_before - val)/smallest_delta*resolution,0)
                    # Fill empty lines
                    for i in range(1, int(delta)):
                        lines.append(empty)
                    lines.append(valform + line)
                value_before = val


            deltas = []
            lines.append([" "*10 + "+" + "-"*50])
        else:
            lines = []

        for i in range(len(lines)):
            print("".join(lines[i]))

    def printCsv(self, reverse=True):
        self.calculateOrbitals()
        lines = []
        for val in self._data:
            for o in self._data[val]:
                line = []
                line.append("%.5f" % (val,))
                # Normalize vector
                o/=np.linalg.norm(o)
                o*=-1 if o[-1] < 0 else 1
                for pos in range(0, len(o)):
                    line.append("%.5f" % (o[-pos-1], ))
                lines.append(','.join(line))

        for i in range(len(lines)):
            if reverse:
                print(lines[i])
            else:
                print(lines[-i-1])




if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='p-Orbitals')

    parser.add_argument('n', metavar='N', type=int,
                   help='Number of orbitals', default=2)
    parser.add_argument('--cyclic', action='store_true', help='Calculate the cyclic variants')
    parser.add_argument('--format', choices=['simple', 'pretty', 'csv'], help='Show eigenvalues and eigenvectors', default='simple')
    parser.add_argument('--reverse', action='store_true', help='Reverses the order of the orbitals (lowest on top)')

    args = parser.parse_args()

    pi = piSystem(args.n, args.cyclic)

    prints = {
        "simple" : pi.printSimple,
        "pretty" : pi.printPretty,
        "csv" : pi.printCsv,
    }

    prints[args.format](args.reverse)