Hückel Molecular Orbital method

huckel.py module

It contains a single class, HMO, implementing the Simple Hückel Molecular Orbital (HMO) method for planar, fully conjugated organic pi-systems.

Given a molecule’s SMILES string, HMO uses RDKit to parse the structure, identify its conjugated pi-system, build the Hückel secular matrix (with Streitwieser or Van Catledge parametrization for heteroatoms, and optional automatic treatment of methyl/alkyl hyperconjugation), and diagonalize it to obtain molecular orbital energies, coefficients, and derived properties (pi charges, bond orders, free valences).

from pyphyschemtools import HMO

benzene = HMO("c1ccccc1")
benzene.summary()
benzene.plot_energy_diagram()

Beyond the core eigenvalue problem, HMO offers a complete toolkit for visualization and analysis: 2D structure display with pi-system highlighting, energy-level diagrams (including a frontier-orbital view), individual or batch molecular-orbital plots on the molecular graph, pi-electron population/charge maps, and export of results to spreadsheets.

Theory

σ/π separability

A planar conjugated molecule is treated in the LCAO-MO framework: each molecular orbital (MO) is written as a linear combination of atomic orbitals (AOs), normally the valence AOs. For a planar system, the AOs split cleanly into two sets with respect to the molecular plane: the in-plane (sp, sp², sp³) AOs building the σ skeleton, symmetric with respect to that plane, and the out-of-plane 2p AOs (one per conjugated atom, perpendicular to the plane) building the π system, antisymmetric with respect to it. Because the two sets transform differently under reflection through the molecular plane, the σ and π problems do not mix: they can be solved separately. Simple Hückel theory only tackles the π part; the σ framework is assumed fixed and is not treated explicitly.

For ethylene, the π AO basis is just the two 2p orbitals perpendicular to the C=C axis, \(\{2p_{y_A}, 2p_{y_B}\}\), combining into a bonding \(\pi\) and an antibonding \(\pi^*\) MO. For benzene, it is the six out-of-plane 2p orbitals, one per ring carbon, \(\{2p_1, \dots, 2p_6\}\): delocalization is simply the π electrons hopping from one such p AO to its neighbors around the ring.

The two Hückel parameters

Solving the π problem means solving the one-electron Schrödinger-like equation \(\hat{h}_\pi \varphi_\pi = \varepsilon_\pi \varphi_\pi\). Rather than evaluate the Hamiltonian matrix elements from first principles, simple Hückel theory introduces two empirical, phenomenological parameters per conjugated system:

• the Coulomb integral \(\alpha_r = \langle 2p_r | \hat{h} | 2p_r \rangle\), the diagonal matrix element of the effective one-electron Hamiltonian for the 2p AO of atom \(r\): in this model it directly is that atom’s orbital energy;

• the resonance (or hopping) integral \(\beta_{rs} = \langle 2p_r | \hat{h} | 2p_s \rangle\): the energy associated with an electron delocalizing between two neighboring 2p AOs \(r\) and \(s\). It is what lets a single secular matrix describe delocalization at all, rather than independent, localized p electrons on each atom.

On top of this, Hückel theory makes two simplifying approximations that turn the general secular problem \(\mathbf{Hc} = E\mathbf{Sc}\) into an ordinary eigenvalue problem:

1. Zero differential overlap: the AO overlap matrix \(\mathbf{S}\) is replaced by the identity, \(\langle 2p_r | 2p_s \rangle = \delta_{rs}\), even though real, neighboring p orbitals do overlap physically;

2. Nearest-neighbor, topological β: \(\beta_{rs}\) is taken to be a single constant \(\beta\) (identical for every formal pi bond) when atoms \(r\) and \(s\) are directly bonded, and exactly zero otherwise, regardless of the actual through-space distance or geometry. The resonance matrix is therefore purely topological: it only encodes the molecule’s connectivity (its sigma-bond graph), not its 3D shape.

Secular matrix and the benzene example

With these two approximations, the π Hamiltonian matrix elements for an \(N\)-atom π system reduce to

\[\begin{split} H_{rr} = \alpha_r \qquad H_{rs} = \begin{cases} \beta & \text{if } r, s \text{ directly pi-bonded} \\ 0 & \text{otherwise} \end{cases} \end{split}\]

(\(\alpha_r = \alpha\) for a plain carbon, shifted for a heteroatom; see below), and the secular equation \(\mathbf{Hc} = E\mathbf{c}\) (\(\mathbf{S} = \mathbf{I}\)) becomes an ordinary matrix eigenvalue problem. For benzene’s six 2p AOs this is, explicitly,

\[\begin{split} \mathbf{H} = \begin{pmatrix} \alpha & \beta & 0 & 0 & 0 & \beta \\ \beta & \alpha & \beta & 0 & 0 & 0 \\ 0 & \beta & \alpha & \beta & 0 & 0 \\ 0 & 0 & \beta & \alpha & \beta & 0 \\ 0 & 0 & 0 & \beta & \alpha & \beta \\ \beta & 0 & 0 & 0 & \beta & \alpha \end{pmatrix} \end{split}\]

the connectivity/adjacency matrix of the ring, dressed with \(\alpha\) on the diagonal and \(\beta\) on each nearest-neighbor off-diagonal entry. Diagonalizing it directly yields the familiar six π MOs and the textbook \(x = 2, 1, 1, -1, -1, -2\) energy pattern (with \(E = \alpha + x\beta\), \(\beta < 0\)).

In practice, the whole problem is made dimensionless by substituting \(x = (E - \alpha)/\beta\), so that diagonalizing \(\mathbf{H}\) amounts to diagonalizing the purely topological matrix \(\mathbf{M}\) (\(\mathbf{H} = \alpha\mathbf{I} + \beta\mathbf{M}\), so \(\mathbf{M}\) is exactly \(\mathbf{H}\) above with \(\alpha \to 0\) and \(\beta \to 1\)): this is exactly what HMO.huckel_matrix() returns, and what the module reports every energy as, in units of \(|\beta|\) relative to \(\alpha\).

Pi population and bond orders

Each occupied MO \(i\), with coefficients \(c_{ri}\) on atom \(r\) and occupation number \(n_i\) (0, 1 or 2), contributes \(n_i c_{ri}^2\) electrons to atom \(r\). Summing over all occupied MOs gives the π electron density on atom \(r\), \(q_r = \sum_i n_i c_{ri}^2\), and comparing it to the number of π electrons that atom formally contributed to the system gives its net π charge. The same coefficients, taken pairwise, give the π bond order between two bonded atoms \(r\) and \(s\), \(p_{rs} = \sum_i n_i c_{ri} c_{si}\), a measure of π bonding character on that bond (close to 1 for a full double bond, close to 0 for a purely single one). Both are standard Hückel population-analysis quantities, and both are available directly from HMO (pi_charge, bond_order, density_matrix(), plot_pi_population()).

Delocalization and conjugation energies

Two complementary measures of pi delocalization are available directly from HMO.

atomization_energy() (and the stored .E_atomization) compares the total pi energy \(E_\pi\) to a reference where every pi electron stays localized on its own atom, non-interacting, at that atom’s own Coulomb level \(\alpha_r\):

\[ E_{at} = E_\pi - \sum_r n_r \alpha_r \]

where \(n_r\) is the number of pi electrons atom \(r\) contributes. For an all-carbon system (\(h=0\) everywhere) this reference is zero, so \(E_{at} = E_\pi\); the gap between the two widens as soon as a heteroatom (\(h \neq 0\)) is involved.

conjugation_energy() (and .E_conjugation) uses a different, more specific reference: a localized Kekulé structure, obtained from RDKit’s own Kekulization of the input SMILES. Each formal pi double bond of that structure is treated as its own isolated, localized 2-center bond: diagonalizing just that \(2\times 2\) block,

\[\begin{split} \begin{pmatrix} \alpha_r & \beta_{rs} \\ \beta_{rs} & \alpha_s \end{pmatrix}, \end{split}\]

and occupying the bonding level with 2 electrons (this reduces to the familiar \(\alpha+\beta\) level for a plain C=C, and generalizes correctly to a heteroatom-involving double bond such as pyridine’s C=N, whose two atoms have different \(\alpha\)). Every atom left unpaired by that Kekulé structure (a lone pair localized on a single center, e.g. pyrrole-type N, furan-type O, a conjugated halogen) keeps the same isolated-atom reference used in atomization_energy():

\[ E_{conj} = E_\pi - \left[ \sum_{\text{double bonds}} 2(\alpha + x_{loc}\beta) + \sum_{\text{unpaired } r} n_r \alpha_r \right] \]

This is the classic textbook “resonance energy” (\(2|\beta|\) for benzene). It is a more specific, usually smaller quantity than the atomization energy, since it isolates the stabilization gained by delocalizing beyond simple, already-bonded, localized double bonds, rather than relative to fully isolated atoms. Cyclobutadiene is the instructive edge case: its localized Kekulé structure (two isolated C=C double bonds) is already as stable as the delocalized, antiaromatic ring, so its conjugation energy comes out to (near) zero. For a non-alternant or unsymmetrical system admitting several valid Kekulé structures, the result follows whichever one RDKit happens to produce, a known limitation of the classical concept itself, not of this implementation.

From the toy model to HMO

Classic Hückel theory, as sketched above, only has one kind of atom (carbon) and one kind of bond. HMO extends this to heteroatom-containing systems the standard way: a heteroatom \(X\) keeps the same secular-matrix machinery, but its Coulomb integral is shifted relative to carbon’s, \(\alpha_X = \alpha_C + h_X \beta_{CC}\), and a C-X bond’s resonance integral is scaled, \(\beta_{CX} = k_{CX} \beta_{CC}\), with the dimensionless \(h_X\)/\(k_{CX}\) pulled from tabulated parameter sets chosen from each heteroatom’s local bonding pattern (pyridine- vs. pyrrole-type N, ether vs. carbonyl O…). Everything downstream (diagonalization, Aufbau/Hund filling, energies, coefficients, charges, bond orders) follows exactly the same recipe as plain benzene above.

Two parameter sets are bundled with the module:

• Streitwieser: A. Streitwieser, Molecular Orbital Theory for Organic Chemists, Wiley, 1961.

• Van Catledge: F. A. Van-Catledge, A Pariser-Parr-Pople-based set of Hückel molecular orbital parameters, J. Org. Chem. 1980, 45 (24), 4801-4802.

API reference

huckel.py

Simple Hückel (HMO) method applied to molecules described by a SMILES string, with Streitwieser parametrization for heteroatoms.

Overview

1. Start from a complete SMILES (sigma skeleton + pi system).

2. RDKit identifies the atoms belonging to the conjugated pi system (aromatic atoms, plus atoms linked by non-aromatic conjugated bonds: dienes, enones, etc.).

3. Each pi atom is given a Coulomb parameter h_X and each pi bond a resonance parameter k_C-X (Streitwieser), relative to the reference values alpha_C and beta_CC:

   alpha_X = alpha_C + h_X * beta_CC
   beta_CX = k_CX   * beta_CC
   

   Automatic detection (degree, number of H atoms, aromaticity) picks the heteroatom type (pyridine-like/pyrrole-like N, ether/furan or carbonyl O, etc.); the caller can always override these values with h_overrides / k_overrides / electron_overrides.

4. The Hückel matrix is diagonalized, the MOs are filled by the Aufbau principle (+ Hund’s rule for degenerate levels), and a number of diagrams / descriptors can be derived from the result: the energy diagram, MO coefficients drawn on the real 2D geometry of the molecule, frontier-orbital (HOMO/LUMO) figures, pi charges and pi bond orders, the density matrix, and reactivity descriptors (HOMO/LUMO energies, hardness, electronegativity…).

Limitations: pure pi-electron model, no explicit self-consistency (no extended Hückel, no charge-dependent VSIP); formal charge only enters through its effect on the number of pi electrons contributed by each atom.

class pyphyschemtools.huckel.HMO(smiles: str, h_overrides: dict[int, float] | None = None, k_overrides: dict[int, float] | None = None, electron_overrides: dict[int, int] | None = None, charge: int = 0, parameter_set: str = 'van_catledge', include_hyperconjugation: bool = True, alkyl_scope: str = 'all', include_halogen_lone_pairs: bool = True)

Bases: object

Hückel pi system built from a SMILES string.

Parameters:

• smiles (str) – SMILES of the complete molecule (sigma skeleton + pi system).

• h_overrides (dict[int, float], optional) – Force h_X for the given RDKit atom index, or k_C-X for the pi bond(s) touching that atom (key = index of one of the two atoms of the bond; use a (i, j) tuple key to target one specific bond).

• k_overrides (dict[int, float], optional) – Force h_X for the given RDKit atom index, or k_C-X for the pi bond(s) touching that atom (key = index of one of the two atoms of the bond; use a (i, j) tuple key to target one specific bond).

• electron_overrides (dict[int, int], optional) – Force the number of pi electrons contributed by the given atom index.

• charge (int, optional) – Defaults to 0. A global formal charge added to the WHOLE pi system, on top of whatever individual atoms already contribute – the easy way to study a system-wide reduction/oxidation (e.g. the benzene dication, HMO(“c1ccccc1”, charge=2), or the naphthalene dianion, HMO(“c1ccc2ccccc2c1”, charge=-2)) WITHOUT having to spell out a carbocation/carbanion center anywhere in the SMILES (which – see _detect_pi_system – also ties that extra/missing electron to one particular atom, rather than letting it delocalize over the whole system the way a real, symmetric di-ion does). Ordinary chemistry sign convention: a positive charge REMOVES that many pi electrons (cation), a negative one ADDS that many (anion); pi_charge’s per-atom values still sum to exactly this charge however it ends up delocalized.

• parameter_set (str, optional) – “van_catledge” (default) or “streitwieser” – which one-/two- center parameter table to use for heteroatoms (see parameter_table). “van_catledge” additionally covers B, Si and P, defines a full heteroatom-heteroatom k_XY matrix (Streitwieser only tabulates k_C-X), and enables free_valence.

• include_hyperconjugation (bool, optional) – Defaults to True. A methyl or alkyl group attached directly to the conjugated system (e.g. toluene’s -CH3, or a benzylic -CH2- of indane/tetralin) is included as one extra pseudo-heteroatom contributing a lone pair (Matsen’s and Stevenson’s “heteroatom model” of hyperconjugation, (sym, kind)=(“C”, “methyl_donor”) – see STREITWIESER_TABLE / VAN_CATLEDGE_TABLE), the whole substituent collapsed to its FIRST attached carbon only (its own further substitution, e.g. the rest of an ethyl or isopropyl group, is not modeled any further). Only a substituent with exactly ONE attachment point to the conjugated system qualifies – a -CH2- bridging two separate pi systems is left alone, as a genuine break in conjugation, not a substituent. The literature value (h=2.0, k=0.7) is strictly only derived for a true methyl; for any larger alkyl it is reused as the closest available approximation, and a warning is raised naming the atom (see alkyl_scope to restrict this to true methyls only, or pass h_overrides/k_overrides/electron_overrides for that atom index to supply a better value by hand). Set to False to restore the previous behavior (these substituents are plain, inert sigma groups, outside the pi system entirely).

• alkyl_scope (str, optional) – “all” (default) or “methyl_only” – with include_hyperconjugation, whether a larger alkyl substituent (-CH2R, -CHR2, -CR3) is included too (reusing the methyl parameter, with a warning) or only a true methyl (-CH3) is. Ignored if include_hyperconjugation=False.

• include_halogen_lone_pairs (bool, optional) – Defaults to True. A halogen (F, Cl or Br) attached directly to the conjugated system (e.g. fluorobenzene’s -F, vinyl chloride’s -Cl) is included in the pi system as one extra 2-electron lone-pair donor, parametrized exactly like any other heteroatom (see STREITWIESER_TABLE / VAN_CATLEDGE_TABLE, (symbol, “halogen”)) – the standard Hückel textbook treatment of aryl/vinyl halides. This is genuinely RDKit’s own blind spot (its GetIsConjugated() is also False for an aryl/vinyl- halogen bond, the same class of issue as the ethylene/allyl cases documented in _detect_pi_system), yet whether to model a halogen’s lone pair at all is a modelling CHOICE rather than an unambiguous correction – some Hückel treatments leave halogens out entirely. A UserWarning is raised for every halogen folded in this way, naming the atom. Set to False to restore the previous behavior (halogens left out, as plain, inert sigma substituents, outside the pi system entirely).

atomization_energy(alpha: float = 0.0, beta: float = 1.0) → float

Hückel atomization (delocalization) energy:

Eat = E_pi(Hückel) - sum_r (alpha_r * contributed_r)

i.e. the total pi-electron energy of the real, delocalized molecule minus the reference energy of the same pi electrons sitting, non-interacting, on their own atoms (each atom r contributing contributed_r electrons at its own Coulomb level alpha_r = alpha + h_r * beta).

Same dimensionless “units of |beta|” convention as total_pi_energy by default (alpha=0, beta=+1, i.e. beta’s magnitude): a positive Eat means the delocalized pi system is more stable than the isolated-atom reference (the usual case: conjugation/aromatic stabilization). Pass the real, physically negative beta (e.g. in eV) for a true physical energy – the sign then flips, since a system that is MORE stable has a LOWER (more negative) energy once beta’s true sign is used.

input: - alpha = Coulomb reference energy (default: 0.0) - beta = resonance reference energy; magnitude convention by default – pass the real, physically negative value for a true physical energy (default: 1.0)

output: - atomization energy Eat, in the chosen alpha/beta units

bond_order(r: int, s: int) → float

Pi bond order between local pi-atom indices r and s (same indexing as self.labels / self.c_MOs rows):

p_rs = sum_i n_i * c_ri * c_si

input: - r = local pi-atom index (row of self.c_MOs) - s = local pi-atom index (row of self.c_MOs)

output: - the pi bond order p_rs between atoms r and s

conjugation_energy(alpha: float = 0.0, beta: float = 1.0) → float

Hückel conjugation (resonance/delocalization) energy, relative to a LOCALIZED (Kekule) reference structure rather than to fully isolated atoms:

Econj = E_pi(Hückel) - [ sum over formal pi double bonds of
                            2*(localized bonding level)
                          + sum over unpaired atoms of
                            (contributed_r * alpha_r) ]

This is the classic textbook quantity (e.g. “benzene’s resonance energy is 2|beta|”): each formal double bond of a Kekule structure (from a Kekulized copy of the input SMILES, see _kekule_pairing) is treated as its own isolated, localized 2-center pi bond – diagonalizing just that 2x2 block (which reduces to the familiar alpha+beta level for a plain C=C, but generalizes correctly to a heteroatom-involving double bond such as pyridine’s C=N, whose two atoms have different alpha) – and occupied with 2 electrons; an atom left unpaired by that Kekule structure (a lone pair localized on a single center: pyrrole- type N, furan-type O, a conjugated halogen, a hyperconjugating methyl/alkyl pseudo-atom, or a formally charged/radical center) keeps the same isolated-atom reference as atomization_energy.

Unlike atomization_energy (reference: fully isolated, non- interacting atoms, no bonding of any kind), this isolates only the stabilization gained by DELOCALIZING beyond simple, already- bonded, localized double bonds – the smaller, more specific number usually meant by “conjugation energy” or “resonance energy” in textbooks.

Same dimensionless “units of |beta|” convention as atomization_energy by default (alpha=0, beta=+1): a positive Econj means the delocalized pi system is more stable than its own localized Kekule structure.

input: - alpha = Coulomb reference energy (default: 0.0) - beta = resonance reference energy; magnitude convention by default – pass the real, physically negative value for a true physical energy (default: 1.0)

output: - conjugation energy Econj, in the chosen alpha/beta units

density_matrix(as_dataframe: bool = True)

First-order pi density (charge/bond-order) matrix:

P = C . diag(n) . C^T

where C holds the MO coefficients (columns) and n the MO occupations. P_rr is the electron density on atom r, P_rs (r != s) is the pi bond order between atoms r and s.

Displayed immediately (works regardless of where this call sits in a cell – see show() for why) and cached as self.density_matrix_df; returns None (read the attribute, or call again, for the value itself).

input: - as_dataframe = if True (default), cache/display the matrix as a labeled pandas DataFrame; if False, as a plain numpy array

output: - None (the matrix is displayed and cached as self.density_matrix_df)

property electron_density: ndarray

Pi electron density per atom (local index, same order as self.labels): q_r = sum_i n_i * c_ri^2.

property free_valence: ndarray

Free valence (radical-reactivity) index per atom:

F_r = Fx_r - sum_s p_rs   (sum over the pi-bonded neighbours s)

where Fx_r is Van Catledge’s reference maximum total bond order for that atom’s type/hybridization (e.g. Fx=1.732=sqrt(3) for an ordinary sp2 carbon). A larger F_r means a more reactive site towards radical attack.

Only defined for parameter_set=”van_catledge” – the Streitwieser table does not provide reference Fx values. Reference: Van- Catledge, J. Org. Chem. 1980, 45 (24), 4801-4802 (Table I).

output: - free valence F_r per atom, as a numpy array (one value per pi atom, same order as self.labels)

huckel_matrix(symbolic: bool = False, as_dataframe: bool = True)

The Hückel matrix itself, such that H = alpha*I + beta*M.

symbolic=False (default): returns the numeric, dimensionless matrix M (diagonal = Coulomb parameters h_r, off-diagonal = resonance parameters k_rs for bonded pi atoms, 0 otherwise). For a simple all-carbon system (all h_r=0, all k_rs=1) this is just 0/1 – e.g. ethylene C=C gives [[0, 1], [1, 0]].

symbolic=True: returns the same matrix with literal Greek-letter alpha/beta entries instead of numbers – a diagonal entry is alpha when h_r=0, else alpha + h*beta with the numeric h_r substituted in; an off-diagonal bonded entry is beta when k_rs=1, else k*beta with the numeric k_rs substituted in; unbonded pairs are 0. Uses real sympy symbols (nicer, typeset rendering in Jupyter) if sympy is installed, otherwise plain Greek-letter strings (“α”, “α + 0.5·β”, …).

Displayed immediately (works regardless of where this call sits in a cell – see show() for why) and cached as self.huckel_matrix_df – a labeled pandas DataFrame (rows/cols = self.labels) if as_dataframe=True and pandas is available, otherwise a raw numpy array (symbolic=False) or a list of lists (symbolic=True). Returns None (read the attribute, or call again, for the value itself).

mo_table(as_dataframe: bool = True, tol: float = 1e-06, with_coefficients: bool = True, coeff_tol: float = 1e-09)

Table of all pi MOs (index, x in units of beta, occupation, role, and – by default – the MO coefficient on every pi atom), ordered from the least stable (top) to the most stable (bottom) MO – same order as the energy diagram. A degenerate HOMO or LUMO (several MOs sharing the same x) is fully reflected in the role column, not just a single representative MO. Values of x within tol of zero are cleaned up to exactly 0.0 (floating- point noise from the diagonalization), and likewise for the coefficients within coeff_tol.

with_coefficients=True (default) adds one column per pi atom (named after self.labels, e.g. “C0”, “N3”, …) holding c_ri, the coefficient of that atom in MO i – set it to False for a compact table with just index/energy/occupation/role (this is what summary() uses, to keep the printed report short for large systems).

Displayed immediately (works regardless of where this call sits in a cell – see show() for why) and cached as self.mo_table_df – a pandas DataFrame if as_dataframe=True and pandas is available, otherwise a list of dicts. Returns None (read the attribute, or call again, for the value itself).

property pi_charge: ndarray

Net formal pi charge per atom = (pi electrons a NEUTRAL atom of that same type would contribute) - (electron density actually found there). Positive means electron-poor relative to the neutral reference, negative means electron-rich.

For almost every atom type this neutral reference IS contributed_pi_electrons (a neutral pyrrole-type N really does contribute 2, a neutral pyridine-type N really does contribute 1 – the heteroatom’s charge state, if any, is already encoded through h instead, e.g. (“N”,”pyridinium”) keeps n_pi=1, same as neutral pyridine). The one exception is a carbon folded into the pi system as (“C”,”anion”)/(“C”,”cation”) (e.g. the allyl/benzyl cation or anion, cyclopentadienyl anion, tropylium): contributed_pi_electrons is 2 or 0 there – the ion’s OWN electron count, needed elsewhere for n_pi_electrons – but the neutral reference used HERE is 1, same as any other carbon, so the molecule’s one formal charge unit is correctly spread over every position by the delocalization itself, instead of being double-counted (pre-paid once at that single atom’s reference, then distributed again by its density). This is why, e.g., tropylium comes out with +1/7 on EVERY carbon (by symmetry), never a spuriously negative nominal “+” center.

Sums to exactly the pi system’s total formal charge (0 for a neutral molecule; +1 for a cation such as tropylium or the allyl cation; -1 for an anion such as the cyclopentadienyl anion) – not always 0.

plot_all_mo(ncols=3, figsize_per=None, bond_length_px=70, save_path=None, dpi=300)

Grid with every pi MO drawn on the molecular graph. See plot_mo_list for a grid restricted to a chosen subset of MOs (e.g. just the frontier ones, or a hand-picked list).

plot_energy_diagram(mode=None, ax=None, figsize=None, tol=1e-06, frontier_figsize=None, save_path=None, dpi=300)

Energy-level diagram. The x axis is oriented so that the energy E = alpha + x*beta increases upward (beta < 0: the largest, most stabilizing x values end up at the bottom).

mode=”FMO” (or “fmo”): instead of the diagram alone, build and return a full figure with the diagram plus the HOMO(s)/LUMO(s) (frontier orbitals) next to it, so the subplots do not have to be composed by hand. Handles degenerate HOMO/LUMO sets and radical (SOMO) cases.

figsize defaults to None, which auto-scales the diagram’s width with the largest level degeneracy (e.g. benzene, max degeneracy 2, stays as compact as before; a highly degenerate spectrum such as C60’s, up to degeneracy 9, gets enough room for every individual segment) – pass an explicit figsize to override.

save_path: if given, the figure is also saved there (see save_fig); dpi is passed through.

input: - mo_index = index of the MO to draw (column of self.c_MOs) - ax = matplotlib Axes to draw on; a new figure is created if None (default) - scale = marker-size scale factor for the largest coefficient - figsize = (width, height) in inches; auto-computed from bond_length_px if None (default) and ax is None - bond_length_px = target on-screen bond length in pixels, used to auto-size figsize (default: 70) - lobe_offset = separation of the two lobes, as a fraction of each atom’s marker radius (default: 0.2); see above - save_path = file path to save the figure to (passed to save_fig); not saved if None (default) - dpi = resolution passed to save_fig when save_path is given (default: 300)

output: - None (the figure is drawn and shown, and saved if save_path is given)

plot_mo_list(mo_indices, ncols=3, figsize_per=None, bond_length_px=70, save_path=None, dpi=300)

Grid with a chosen list of pi MOs drawn on the molecular geometry (see plot_mo_on_graph for a single one) – e.g. mol.plot_mo_list([0, 1, 7, 8]) to compare a few MOs of interest side by side, in exactly the order given, without plotting the whole set as plot_all_mo does.

figsize_per (the per-panel size) defaults to None, which auto-sizes EVERY panel, uniformly, from the molecule’s own 2D layout so a C-C bond renders at roughly bond_length_px screen pixels – same idea, and the same _figsize_for_bond_px helper, as the standalone plot_mo_on_graph’s own figsize/ bond_length_px (all panels here show the SAME molecule, just a different MO each, so one computed size applies to all of them). Pass an explicit figsize_per to override.

save_path: if given, the figure is also saved there (see save_fig); dpi is passed through.

input: - mo_index = index of the MO to draw (column of self.c_MOs) - ax = matplotlib Axes to draw on; a new figure is created if None (default) - scale = marker-size scale factor for the largest coefficient - figsize = (width, height) in inches; auto-computed from bond_length_px if None (default) and ax is None - bond_length_px = target on-screen bond length in pixels, used to auto-size figsize (default: 70) - lobe_offset = separation of the two lobes, as a fraction of each atom’s marker radius (default: 0.2); see above - save_path = file path to save the figure to (passed to save_fig); not saved if None (default) - dpi = resolution passed to save_fig when save_path is given (default: 300)

output: - None (the figure is drawn and shown, and saved if save_path is given)

plot_mo_on_graph(mo_index, ax=None, scale=800, figsize=None, bond_length_px=70, lobe_offset=0.2, save_path=None, dpi=300)

Coefficients of MO mo_index drawn on the real 2D geometry of the molecule, each pi atom rendered as a stylized two-lobed p orbital rather than a flat disc: a front lobe (circle size = |c_i|, color = sign of c_i, same convention as before) shifted slightly toward the bottom-left along the diagonal, and a same- size, opposite-color back lobe shifted the same amount toward the top-right, drawn behind it. The atom label stays at the true atom position, between the two lobes.

figsize defaults to None, which auto-sizes the panel from the molecule’s own 2D layout so a C-C bond renders at roughly bond_length_px screen pixels (see _figsize_for_bond_px) – this is what keeps ethylene’s one short bond from filling the whole panel the way it would at a fixed figsize, while a much larger pi system (coronene…) still gets the room it needs; pass an explicit figsize to override. Only applies when this call creates its OWN figure (ax=None) – when composed into an existing ax (plot_all_mo/plot_mo_list/the frontier- orbital panels), the panel’s size is already decided by the caller instead.

lobe_offset: how far apart the two lobes are drawn, as a fraction of EACH atom’s own rendered circle radius (so a big coefficient’s lobes are offset more than a small one’s, and the two always overlap sensibly regardless of marker size) – 0 collapses back to a single disc per atom (the pre-existing look); 1 offsets each lobe by a full radius (touching at the atom’s true position); the default (0.2) gives a subtle, mostly-overlapping two-lobed look – raise it for a more pronounced dumbbell shape.

save_path: if given, the figure is also saved there (see save_fig); dpi is passed through.

input: - mo_index = index of the MO to draw (column of self.c_MOs) - ax = matplotlib Axes to draw on; a new figure is created if None (default) - scale = marker-size scale factor for the largest coefficient - figsize = (width, height) in inches; auto-computed from bond_length_px if None (default) and ax is None - bond_length_px = target on-screen bond length in pixels, used to auto-size figsize (default: 70) - lobe_offset = separation of the two lobes, as a fraction of each atom’s marker radius (default: 0.2); see above - save_path = file path to save the figure to (passed to save_fig); not saved if None (default) - dpi = resolution passed to save_fig when save_path is given (default: 300)

output: - None (the figure is drawn and shown, and saved if save_path is given)

plot_pi_population(ax=None, figsize=None, bond_length_px=100, charge_fmt='{:+.2f}', bond_fmt='{:.2f}', atom_fontsize=10, bond_fontsize=10, save_path=None, dpi=300)

Classic Hückel ‘charge and bond-order diagram’: the pi system graph annotated with the net pi charge on each atom and the pi bond order on each bond (see pi_charge / bond_order / density_matrix).

figsize defaults to None, which auto-sizes the figure from the molecule’s own 2D layout so a C-C bond renders at roughly bond_length_px screen pixels (see _figsize_for_bond_px) – larger than plot_mo_on_graph’s own default (100 vs. 70) since the bond-order labels sitting right on top of each bond need more room to stay legible than a plain p-orbital lobe does; pass an explicit figsize to override. atom_fontsize/ bond_fontsize control the label sizes directly – bump them further for a very crowded graph, or when exporting to a larger figure.

save_path: if given, the figure is also saved there (see save_fig); dpi is passed through.

input: - mo_index = index of the MO to draw (column of self.c_MOs) - ax = matplotlib Axes to draw on; a new figure is created if None (default) - scale = marker-size scale factor for the largest coefficient - figsize = (width, height) in inches; auto-computed from bond_length_px if None (default) and ax is None - bond_length_px = target on-screen bond length in pixels, used to auto-size figsize (default: 70) - lobe_offset = separation of the two lobes, as a fraction of each atom’s marker radius (default: 0.2); see above - save_path = file path to save the figure to (passed to save_fig); not saved if None (default) - dpi = resolution passed to save_fig when save_path is given (default: 300)

output: - None (the figure is drawn and shown, and saved if save_path is given)

properties(alpha: float = 0.0, beta: float = 1.0, as_series: bool = False)

Frontier-orbital energies and simple conceptual-DFT reactivity descriptors, handy to collect into a table (DataFrame) across several molecules.

alpha/beta control E_HOMO, E_LUMO, gap, hardness, softness, chemical_potential, electronegativity and electrophilicity. They default to alpha=0, beta=+1: E_HOMO, E_LUMO, chemical_potential and electronegativity are then reported as their COEFFICIENT of beta (E = alpha + value*beta), the same convention as mo_table (x_beta) – e.g. E_HOMO = 0.618 means alpha + 0.618*beta. gap, hardness and electrophilicity are intrinsically positive quantities and are always returned as positive values, in units of |beta| (softness in 1/|beta|). Pass the real, physically negative beta (e.g. in eV) to get all of these in physical units.

total_pi_energy, atomization_energy and conjugation_energy are NOT affected by alpha/beta here – they are always reported in the dimensionless “units of |beta|” convention (beta’s magnitude, always positive for a normal bonded pi system; see total_pi_energy/atomization_energy/conjugation_energy, or self.E_total / self.E_atomization / self.E_conjugation). Call those methods directly, with your own alpha/beta, for a physical-unit total/atomization/conjugation energy.

input: - alpha = Coulomb reference energy (default: 0.0) - beta = resonance reference energy; dimensionless magnitude convention by default – pass the real, physically negative value (e.g. in eV) for physical units (default: 1.0) - as_series = if True, return a pandas Series instead of a dict (requires pandas); default: False

output: - dict (or pandas Series if as_series=True) with: smiles, n_pi_atoms, n_pi_electrons, homo_index, lumo_index, E_HOMO, E_LUMO, gap, hardness, softness, chemical_potential, electronegativity, electrophilicity, total_pi_energy, atomization_energy, conjugation_energy

save_results(path_with_name: str, tol: float = 1e-06, **kwargs)

Save the MO table (with coefficients) and the density matrix together in one multi-sheet Excel/ODS file (one sheet each: “MO”, “density_matrix”), via pandas.ExcelWriter – the engine (openpyxl for .xlsx, odfpy for .ods) is picked automatically from the file extension, and must be installed. For a single table (e.g. properties_table’s own result), just call the module-level save_data(df, path) directly instead.

tol: passed to the MO table (numerical cleanup of near-zero values – not the coefficient tolerance, which keeps its own default; pass coeff_tol via a direct mo_table() call if you need to change that one too).

Computes the two tables directly (rather than going through mo_table/density_matrix) so saving a file stays silent – those public methods now always display their result immediately, which isn’t wanted here. Both tables are still cached as self.mo_table_df/self.density_matrix_df, same as calling each method would.

show(ax=None, figsize=None, size=None, bond_length_px=80, title=None, show_indices=True, show_hybridization=False, highlight_pi=True, save_path=None, dpi=300)

Draw the molecule’s FULL 2D skeleton (sigma + pi, with implicit H’s, wedge/dash stereo if any…) using RDKit’s own depiction engine – unlike plot_mo_on_graph / plot_pi_population, which only draw the pi subsystem detected for the Hückel treatment. Handy as a quick sanity check that the SMILES was parsed/drawn (and the conjugated system detected) as intended.

show_indices (default True): label each atom with its RDKit atom index – the same index used in self.pi_atom_idx and in the “C0”, “N3”… entries of self.labels. show_hybridization: annotate each atom with its hybridization (sp, sp2, sp3…) as given by RDKit. highlight_pi (default True): highlight the atoms/bonds of the conjugated pi system actually used for the Hückel treatment (self.pi_atom_idx) – the rest of the molecule (sigma-only substituents, sp3 carbons…) is left unhighlighted.

size defaults to None, which sizes RDKit’s own drawing canvas from the molecule’s actual 2D layout (its full sigma+pi skeleton, not just the pi system) so a C-C bond renders at a FIXED bond_length_px (80 by default) – RDKit’s own fixedBondLength drawing option (verified empirically: it is NOT a 1:1 pixel setting on this RDKit build, so the conversion is calibrated against _bond_unit_length). This is what keeps ethylene’s one bond from filling the whole canvas the way a fixed size would, while a large, crowded molecule (coronene…) still gets a canvas sized to actually fit it, hybridization labels and all. Pass an explicit size (pixels, RDKit’s own legacy behavior: the molecule is auto-scaled, in either direction, to fill that fixed canvas) to go back to a plain fixed-size drawing instead – bond_length_px is then ignored. The drawing’s blank margin is cropped automatically either way.

figsize defaults to None, which sizes the matplotlib figure to match the (now-cropped) drawing 1:1 in pixels when size was auto-computed (so the bond-length guarantee above survives all the way to the actual on-screen rendering, not just the intermediate RDKit image), or to the drawing’s aspect ratio at a fixed 4-inch height when an explicit size was given instead (the previous behavior). Pass an explicit figsize to override either way (e.g. to make it larger, or to match a figure you’re composing by hand). title defaults to the SMILES. save_path/dpi: see save_fig.

summary()

Print a short text report: SMILES, pi atoms (with the number of pi electrons each one contributes), total pi electrons, total Hückel pi-energy, atomization energy, conjugation (Kekule- referenced resonance) energy, and the full MO table (see mo_table to get it as a DataFrame instead of printed text).

total_pi_energy(alpha: float = 0.0, beta: float = 1.0) → float

Total Hückel pi-electron energy:

E_pi = sum_i n_i * (alpha + x_i * beta)

where n_i is the occupation of MO i and x_i its coefficient (self.E_MOs). Defaults to the dimensionless “units of |beta|” convention used when quoting these energies in textbooks (e.g. “E_pi(benzene) = 6*alpha + 8*beta”): alpha=0, beta=+1 – the magnitude of beta, reported as a plain coefficient, not beta’s true physical sign. Since the lowest-energy (most bonding) MOs are filled first and have the largest x, this dimensionless total is positive for any normal bonded pi system. Pass the real, physically negative beta (e.g. in eV) to get the true, physical (and then negative) energy.

input: - alpha = Coulomb reference energy (default: 0.0) - beta = resonance reference energy; magnitude convention by default – pass the real, physically negative value for a true physical energy (default: 1.0)

output: - total pi-electron energy E_pi, in the chosen alpha/beta units

used_parameters(as_dataframe: bool = True)

The parameter-table rows actually used for the atom classes present in THIS molecule – a convenient subset of parameter_table(self.parameter_set) to sanity-check what was applied (one row per distinct (symbol, kind), in the order first encountered). Each row also repeats self.parameter_set (so the table is self-identifying – “streitwieser” or “van_catledge” today, whatever other parameter sets may be added later – even when read out of context, e.g. printed on its own or saved to a file).

Displayed immediately (works regardless of where this call sits in a cell – see show() for why) and cached as self.used_parameters_df – a pandas DataFrame if as_dataframe=True and pandas is available, otherwise a list of dicts. Returns None (read the attribute, or call again, for the value itself).

pyphyschemtools.huckel.hmo_properties_table(smiles_list, alpha: float = 0.0, beta: float = 1.0, **huckel_kwargs)

Build a table of frontier-orbital descriptors (see HMO.properties) for a list of molecules, ready to use as a pandas DataFrame.

Parameters:

• smiles_list (list) – Each element is either a bare SMILES string, or a (smiles, overrides_dict) tuple for per-molecule keyword overrides (h_overrides, k_overrides, electron_overrides, charge).

• alpha (float) – Passed to HMO.properties for every molecule.

• beta (float) – Passed to HMO.properties for every molecule.

• **huckel_kwargs – Default keyword arguments passed to every HMO(…) construction (overridden per-molecule by the tuple form above).

• available (Returns a pandas DataFrame if pandas is)

• list (otherwise a)

• dicts. (of)

pyphyschemtools.huckel.parameter_table(parameter_set: str = 'van_catledge', as_dataframe: bool = True)

The full one-center Hückel parameter table for parameter_set (“streitwieser” or “van_catledge”): one row per (symbol, kind) with its h, k_to_carbon and n_pi_electrons, plus Fx (reference free valence) for “van_catledge” (Streitwieser does not define one).

Returns a pandas DataFrame if as_dataframe=True and pandas is available, otherwise a list of dicts.

References

streitwieser: A. Streitwieser, Jr., “Molecular Orbital Theory for

Organic Chemists”, John Wiley & Sons, 1961, p. 135.

van_catledge: F. A. Van-Catledge, “A Pariser-Parr-Pople-Based Set of

Hückel Molecular Orbital Parameters”, J. Org. Chem. 1980, 45 (24), 4801-4802.

pyphyschemtools.huckel.save_data(df, path_with_name: str, **kwargs)

Saves a Pandas DataFrame to CSV or Excel, automatically creating missing directories.

Parameters:

• df (pd.DataFrame) – The DataFrame to save.

• path_with_name (str) – The path and filename (e.g., “data/results.csv”).

• **kwargs – Additional arguments passed to to_csv or to_excel (e.g., index=True).

pyphyschemtools.huckel.save_fig(path_with_name: str, fig=None, dpi: int = 300, **kwargs)

Saves a figure to a path, automatically creating missing directories. Works with the current plt or a specific figure object.

Parameters:

• path_with_name (str) – The path and filename (e.g., “results/plots/energy.svg”).

• fig (matplotlib.figure.Figure, optional) – A specific figure object. If None, uses plt.gcf().

• dpi (int) – Resolution for raster formats. Defaults to 300.

• **kwargs – Additional arguments passed to savefig (e.g., transparent=True).