The orbital layer

Every chemical picture in this workbench comes from one Python package, qo_chem, run on one 3D conformer. The 3D viewer's hybrids, the π density and pKa maps, the charges and bond orders in /api/model, and the forecasts in /api/pi-model are all read from the same result. If two pages disagree, one of them has a bug: there is no second model to blame.

The layer is deliberately small. It uses chemical rules that each have a name and a test behind them, so you can follow every number back to its origin. It is not Hartree–Fock or DFT, and it says so wherever it is shown.

#The pipeline

SMILES
  → RDKit: graph, Kekulé bonds, one ETKDG conformer (MMFF94 or UFF)
  → normalize      hypervalent X=O → X⁺–O⁻, metals cut into ion pairs
  → electronic     which atoms join the π system, σ and π lone pairs
  → huckel         π orbitals with p axes from the geometry, per component
  → radial, basis  Hartree–Fock valence functions, 3D overlaps, Löwdin basis M = S^-1/2
  → hybrids        σ hybrids that follow the real bond angles
  → pka, protonation, acceptors, hammett    chemistry built on the above
  → provenance     SHA-256 of the result, for reproducibility

One call, qo_chem.analysis.analyze(graph), runs all of it and returns the chemistry block of /api/model.

#1. Which electrons are π electrons?

The gut feeling. In a carboxylate both oxygens are the same, whatever Kekulé structure the SMILES happened to draw. Any model that gives them different charges is reading the drawing, not the molecule.

The rule. An atom's role is decided from its σ-bond count and its Lewis lone pairs, never from where the SMILES put the double bond:

A donor oxygen has one lone pair in the p orbital that overlaps the ring (the π lone pair) and one in the bond plane (the σ lone pair). Each π atom gets Hückel parameters from Streitwieser's table (Molecular Orbital Theory for Organic Chemists, 1961): \(\alpha_X = \alpha + h_X\beta\) and \(\beta_{CX} = k_{CX}\beta\), with \(h = 1.0\) for O= and 2.0 for a donor O:, 0.5 for N= and 1.5 for N:.

The test. For all 20 molecules of the pKa study, the automatic typing reproduces the hand-assigned types, and both oxygens of every carboxylate and nitro group come out with identical π charges.

#2. π orbitals on the real geometry

Hückel (1931) assumed a flat molecule with every p orbital parallel. Real molecules twist. Here each π atom gets its own p axis \(\hat n_r\) from the conformer: the normal of its three bonds for a trigonal centre, two perpendicular axes on an sp centre, and the lone-pair direction on a pyramidal donor.

The coupling between two bonded orbitals is scaled by how well their axes line up across the bond,

A nitro group rotated 90° out of the ring therefore decouples by itself; no special case is written for it. A pyramidal donor such as aniline's nitrogen couples with amplitude \(\sqrt{1-s}\) of its lone-pair hybrid: the more s character the lone pair keeps, the less of it is available to the ring. A centre counts as partial when the three bond angles add up to 335–355° (360° is flat, 328° is tetrahedral).

Diagonalising the Hückel matrix of each conjugated component gives orbital energies \(E = \alpha + x\beta\), coefficients, occupations and the density matrix P, whose diagonal holds the π populations and whose off-diagonal holds Coulson's bond orders.

Check: the π bond order of benzene is 2/3, to 1 part in 10⁴. The server repeats this test at start-up and reports it at /api/pi-density-health.

#3. Orbitals you can draw: Hartree–Fock radial functions and Löwdin's basis

To draw a density you need functions in space, not just coefficients. Each π orbital is a valence p function of its atom, pointing along its own axis,

Which radial function R? Until version 16.0 the layer used Slater's rules (1930), one exponent per element. An experiment against Hartree–Fock atoms (Radial functions: an experiment) showed that Slater's 2p is about 12 % too compact and gives π overlaps 26–38 % too small. Since 16.1 every element uses its free-atom Hartree–Fock valence orbital, written as a short sum of Slater functions (1–6 terms, overlap with Hartree–Fock ≥ 0.9999),

the idea of Clementi and Roetti (1974), refitted for the elements the app uses (qo_chem/data/radial_fits.json, generated by tools/radial/make_radial_table.py). The same table drives the 3D viewer's hybrid surfaces and the π maps.

Two p orbitals with arbitrary axes overlap as

where u points from one nucleus to the other. For sums of Slater functions both radial integrals are sums of Slater-pair integrals, done by quadrature in prolate spheroidal coordinates. Checks: the term sums agree with a dense direct integration to 10⁻⁷; for two carbons 1.39 Å apart \(S_\pi = 0.334\) (Slater's rules gave 0.248), and \(S_\sigma = 0.226 < S_\pi\), because the Hartree–Fock lobes reach past each other at that distance.

What this changes: the π overlaps, the Löwdin matrix and therefore the density, interference and orbital maps. What it does not change: Hückel energies, charges, bond orders and the pKa descriptors, which never use a radial function.

Hückel theory pretends the overlap is zero. That is exactly true in Löwdin's symmetric basis (1950), \(M = S^{-1/2}\), so the density in space is

It splits into an atomic part, \(\sum_r P_{rr}\lambda_r^2\) (the λ_r are the Löwdin orbitals), and an interference part, the rest. Because Löwdin orbitals are orthonormal, the interference part integrates to exactly zero: it moves charge into the bonds without creating any. That is Ruedenberg's picture of the chemical bond (1962), and it is what the "interference" view of the π density map shows.

#4. σ hybrids that follow the bond angles

Textbooks say water's oxygen is sp³. Its bond angle says otherwise. Coulson (Valence, 1952) showed that two hybrids \(h = s + \sqrt{\lambda}\,p\) at angle θ are orthogonal only when

so at 104.5° the O–H hybrids are sp\(^{3.99}\), not sp³. The layer builds the hybrids from the conformer's real bond vectors: seeds are made orthonormal with Löwdin's procedure and the s weights are fitted (Levenberg–Marquardt) so each hybrid points as closely as possible along its bond. Where no exact solution exists, as in three- and four-membered rings, the least-bent set comes out and the bend is reported.

Lone pairs can be shown as a σ + π pair (default, it matches photoelectron spectra) or as an equivalent "rabbit-ear" pair; the two are a rotation of each other and describe the same electrons.

#5. Acidity from the π layer

For each ionisable site the acid and its conjugate base are built on the same geometry: only the proton moves. The layer is run on both, and three descriptors are read off:

The forecast lines were fitted on the curated set of /pka-density and are refitted on the server's own conformers by tools/calibrate_pka.py:

An RMSE of about 1.3 pKa units means "ranks substituents correctly and gets the size of the effect", not "replaces a measurement". The amine line rests on five molecules: treat it as a hint.

Around the forecasts sit three well-tested empirical tools:

• Group pKa values with Henderson–Hasselbalch give the dominant form at a chosen pH (protonation). These are reference values, not predictions.
• Hammett constants σ_m and σ_p for 28 substituents (Hammett 1937; values from Hansch, Leo and Taft 1991) give pKa = 4.20 − Σσ for substituted benzoic acids. Ortho substituents are refused: σ does not describe them.
• Acceptor orbitals. σ* of an sp³ C–X bond from a two-centre model with Pauling electronegativities (for CH₃Cl: 0.80 on C, −0.60 on Cl, so the nucleophile attacks carbon from the back), and π* from the Hückel LUMO, with the Bürgi–Dunitz approach angle of 107° on carbonyl carbons (1974).

#6. Molecules the simple rules would get wrong

• Hypervalent atoms. A sulfone or phosphate "S=O" is not a π bond (there is no d-orbital participation). The graph is charge-separated to X⁺–O⁻ before the π layer sees it; the SMILES you typed is not changed.
• Metal salts. Na–O and similar bonds are cut into ion pairs; ions are drawn as empty ns orbitals (Lewis acids) with an optional water shell (Mg–O 2.07 Å).
• Ions without a π system simply have no π component. The layer reports this instead of inventing one.

#7. Reproducibility

Every result carries a provenance record: the inputs (SMILES, seed, pH, force field), the package version, and a SHA-256 hash of the canonical JSON with floats rounded, so noise in the last bit does not change the hash but any chemical change does. The same SMILES and seed on the same version give the same hash.

#8. Cost and scaling

Until 16.4 the ring search ran a breadth-first search from every bond, through the whole molecule, once per atom. That was invisible at 30 atoms and took 71 s on a 384-atom protein fragment. Limited and cached, the same fragment takes 1.1 s, and 5000 atoms take 5.5 s (tools/sites/RESULTS.md). Zinc joined the metals in 16.4: Zn²⁺ is an ion whose empty 4s/4p (virtual Hartree–Fock orbitals of Zn²⁺) accept the ligands' pairs. At a protein site they are drawn as acceptor hybrids (qo_chem/site_orbitals.py, Protein sites).

#9. What this model leaves out

• No self-consistency: charges do not feed back into the orbital energies.
• σ inductive effects enter only through Hammett constants and electronegativity, not through the π charges.
• No solvent: pKa values in water depend on solvation that is not computed here.
• One conformer: flexible molecules have many.
• Heteroatom parameters are empirical; they are good for first-row atoms and rougher for S, P and the heavier halogens (marked as approximate in the output).
• The radial functions are those of free atoms; orbitals contract somewhat in molecules.

When a result depends on one of these, the API adds a warning or a badge rather than hiding it.

#Where to go next

• The maps built on this layer: π density field engine and the page π electron density and pKa maps.
• The JSON it produces: Using the API.
• The code: qo_chem/ in the server package, with its tests in tests/test_qo_chem.py.