How the orbital diagrams are computed

The orbital diagrams of CrystOD — the molecular ones of crystod-mol --diagram (section 33) and the crystal ones of crystod --diagram (section 3) — are built by the same recipe: symmetry-adapt the fragment orbitals, evaluate every overlap between them, and solve the resulting eigenvalue problem. This page describes that engine, in its qualitative (extended-Hückel) and quantitative (PySCF) forms. Reading the diagrams themselves needs none of this; the command pages are self-contained.

How the molecular-orbital diagram is built (extended-Hückel engine)

The construction follows the textbook route, made quantitative step by step:

1. the molecular point group is detected and the ligand atomic orbitals are symmetry-adapted per irrep (the molecular SALCs of crystod-mol --element --orbital, section 32) — for NH3, A1: [1s(H1) + 1s(H2) + 1s(H3)] and E: [1s(H1) - 1s(H3), 1s(H1) - 2 1s(H2) + 1s(H3)];

2. every valence orbital — the outer, chemically active orbital of each atom — is represented by a single-zeta Slater-type orbital (STO), i.e. by one exponential function per orbital, with the standard extended-Hückel exponents. Each orbital also carries a diagonal energy H_ii, taken from its valence-state ionization energy, which measures how tightly that orbital binds an electron in the free atom;

3. all two-center overlap integrals between the STOs — the numbers S that measure how much two orbitals on different atoms occupy the same region of space — are evaluated exactly, by Gauss-Laguerre x Gauss-Legendre quadrature in prolate-spheroidal coordinates, to machine precision and validated against the analytic H 1s-1s formula. Ligand-ligand overlap is therefore not neglected. The overlaps between each ligand SALC and each central-atom atomic orbital (AO) are printed per irrep (NH3: < a1 | N 2s > S = 0.72, < e | N 2p > S = 0.57);

4. orbitals belonging to different irreps cannot mix, so the problem falls apart into one small block per irrep. Within each block the generalized eigenvalue problem — the secular equation whose solutions are the molecular-orbital energies and their coefficients — is solved with the Wolfsberg-Helmholz off-diagonal elements H_ij = K S_ij (H_ii + H_jj)/2 (K = 1.75). A large overlap integral therefore directly produces a large splitting between the bonding and the antibonding level (this is symmetry-adapted extended Hückel).

Why ghost atoms? The basis-set superposition error (BSSE)

A quantum-chemistry calculation builds the wave function — the mathematical object from which every property of the electrons follows — out of a finite set of atom-centred basis functions, so its quality depends on how many functions sit near the electrons. Compute the O atom of H2O in the basis of O alone, then the whole molecule in the basis of O and the two H atoms, and the molecule has extra freedom to describe its electrons — more functions to build them from — that the isolated fragment never had. Its energy drops for a reason that has nothing to do with chemical bonding. The spurious part of the stabilization is the basis-set superposition error (BSSE) — an artefact of the finite basis rather than a physical effect. It makes computed interaction energies look too attractive, and it does not vanish simply by enlarging the basis on one side only.

The standard cure is the counterpoise correction of Boys and Bernardi: compute each fragment in the full molecular basis, keeping the basis functions of the removed atoms but not their nuclei or electrons. Such nucleus-free, electron-free basis centres are called ghost atoms. Both fragments then carry exactly the same basis-set quality as the molecule, so the difference E(mol) − E(left) − E(right) contains bonding only.

CrystOD needs this twice over. First, the interaction energy printed under the diagram is BSSE-free. Second, because all three calculations live in one and the same AO space, every molecular MO can be projected exactly onto the fragment MOs — written as an exact combination of them, with no fitting and no approximation. That is what makes the correlation lines and the percentage compositions quantitative rather than merely indicative. The price is that a few fragment solutions end up living mostly on the ghost centres; those are filtered out by the real-atom Mulliken-population test described on the crystod-mol --diagram --pyscf command page (section 33).

Two conventions that keep the molecular fragment columns honest

Spin and spherical averaging. The electrons of a partly filled frontier shell are spread evenly over it, so every degenerate orbital of the shell takes the same fractional occupation. Without this, a partially filled shell (C 2p², the t2² of an H4 cage) would break the point-group symmetry of the fragment, whereas the pre-bonding reference is supposed to preserve it.

The ghost-population filter. A fragment level may end up living mostly on the ghost centres — a variational tail of the fragment density toward the atoms that were removed, not a genuine state of that fragment. Each fragment level is therefore judged by its real-atom Mulliken population, the share of the density that Mulliken’s partitioning assigns to the real atoms rather than to the ghost centres. Below 35% the level is a BSSE artefact and is dropped from the diagram and from the compositions; the surviving projections are weighted by the real-atom population. In benzene this is what makes the pure-π 1e1g HOMO come out as 100% C6 1e1g, with no spurious “H6” character borrowed from ghost-carbon levels. Fragment orbital sketches likewise draw the real atoms only, and degenerate partners are canonicalized so that they match the SALCs drawn by --visualize exactly.

In the crystalline --pyscf diagram the same tail appears in the irreps where the fragment has a partner (for ScF3: GM4⁻ = Sc p, R1⁺ = Sc s, and not in GM5⁻); there it is reported as the ghost weight per level rather than filtered, and --no-ghost removes the ghost basis from the fragment variational space entirely — at the cost of the counterpoise consistency (~2 mHa BSSE per ScF3 fragment).

How the crystal-orbital diagram is built (extended-Hückel engine)

This is what crystod --diagram --co-left ... --co-right ... (section 3 of crystod) evaluates at every special k point of the space group:

1. the Bloch orbitals of each fragment (e.g. Sr 1s…5p + Ti 1s…3d | O 1s 2s 2p) are symmetry-adapted per irrep of the little group of k (the site-symmetry induced representations of the --atomic-orbital hybridization analysis);

2. all intra- and inter-fragment overlaps are evaluated exactly as Bloch lattice sums of STO overlap integrals (single-zeta s/p, standard double-zeta d, sigma/pi/delta Slater-Koster assembly rotated with the exact real Wigner-D matrices), with a per-shell-pair lattice-sum cutoff probed from the actual STO tails (diffuse cation shells reach 30+ bohr);

3. the point-charge ligand field enters as exact same-site matrix elements <phi_i|q/|r-R||phi_j> — Laplace expansion into real spherical harmonics (the wigner_D_real real-orbital convention of the representation-matrix background above) with closed-form radial integrals, Loewdin-consistent with the identity on-site metric; near charge shells explicitly, the long-range tail by Ewald summation — validated against the NaCl Madelung constant and the exact 6Dq/-4Dq octahedral splitting ratio, so the fragment d states show the true electrostatic t2g/eg splitting. The background-dependent monopole of the (charged) sublattice array is omitted: it largely cancels against the intra-atomic charging energy absent from the extended-Hueckel VSIPs, and the omitted jellium-referenced values are printed for comparison with charged-cell DFT references;

4. the generalized eigenvalue problem with the Wolfsberg-Helmholz Hamiltonian H_ij = K S_ij (H_ii + H_jj)/2 (the diagonal carries the same-orbital neighbour-cell Bloch sums; the ligand-field blocks are added identically to the fragment and crystal Hamiltonians, so the three columns share one energy reference, and the symmetry invariance of H is self-checked per k point) is solved: fragment orbitals sharing an irrep split into bonding and antibonding crystal orbitals, orbitals without a partner remain rigorously nonbonding — the COD mixing rule.

The interactive HTML page that these steps produce is described with the command in section 3 of crystod.

One numerical caveat belongs to step 4. In dense cation sublattices the diffuse valence shells (Sr 5s/5p, Ti 4s/4p, …) make a few Bloch combinations nearly linearly dependent on the rest of the basis, and for those the extended-Hückel energies diverge — the well-known overlap catastrophe. Combinations with an overlap eigenvalue below 0.2 are therefore removed by canonical orthogonalization, and the terminal report counts them per k point. A second, purely conventional caveat: the irrep labels at zone-boundary points depend on the origin of the input structure, as in every SALC analysis, so compare against a reference in the same setting.

Ghost atoms, cell neutrality and the deep-level alignment (crystal --pyscf)

In the quantitative crystal-orbital diagram (crystod --diagram --pyscf, section 3 of crystod) the three periodic SCF calculations — left fragment, right fragment, crystal — are tied together like this.

The removed sublattice stays in the basis as ghost atoms (basis functions without a nucleus), so all three calculations span the same AO space and the crystal orbitals can be projected exactly onto the fragment Bloch orbitals (counterpoise-consistent, as in the molecular version); in addition it acts on the fragment as its formal-charge point lattice, evaluated as a jellium-referenced FFT potential — validated against crystod.point_charge_field.ewald_site_potential to 2e-6 eV for a neutral charge array. Each fragment is therefore one sublattice in the Madelung field of the other: the electronic state before chemical bond formation.

Both the remaining ions and the replacing point charges carry the formal oxidation states (--oxidation Sc=+3 F=-1 to override), which makes every cell neutral: the monopole divergence disappears, every electron count is even, and the fragment counts add up to the crystal’s (ScF3: 8 + 24 = 32). Neutrality alone does not put the three calculations on one absolute scale, though — each pins its own G = 0 (cell-averaged) potential to zero, and the value of that average depends on the second moment of the cell’s charge density, which changes when an ion is replaced by a bare point charge. The raw columns are therefore offset by one rigid, k-independent constant each (ScF3: Sc column ~2 eV, F3 column ~5 eV below the crystal), the same reference problem as in band-offset calculations. The diagram removes it by deep-level (XPS-style) alignment: the deepest chemically inert fragment level — one that reappears in the crystal at ≥ 80% purity, e.g. the Sc 3s semicore — must keep its energy across bond formation, and the energy zero is put at its pre-bonding (fragment) value. Only the near-purest anchor pairs across the k points set the constant (printed with anchor, purity and k-spread); --no-align keeps the raw references instead.

Why --onsite takes the fragment columns from the crystal Fock

The default crystod --diagram --pyscf mode (section 3 of crystod) solves each fragment as its own SCF (formal-charge ions + point charges + ghost basis), so the three columns live on three different electrostatic references; deep-level alignment removes one rigid constant per column, but the site-dependent part of the point-charge-model error remains — a weakly bonding level can appear ~1.5 eV above its parent (SrTiO3: the chemically inert Sr 4s column proves the shift is pure environment). --onsite removes that error structurally, by reading the fragment columns off the converged crystal Fock operator instead.

Whole-sublattice blocks are not used for those columns: with a diffuse basis the raw AO block lets cation functions fall variationally into the removed side’s potential wells (disguised anion states), and orthogonalized blocks load strongly-overlapped shells with large orthogonalization penalties — the per-shell Rayleigh quotient has no variational freedom to abuse.

What the composition percentages mean (Löwdin vs Mulliken)

The composition list quoted in three places — the Level-details bars, the hover tooltip and the terminal — is one and the same list: the per-(element, shell) AO population of the eigenvector, the same partial-charge measure as crystod --dos --pyscf. The default is Löwdin, |c|² in the symmetrically orthogonalized basis: the orthonormal set closest to the atomic orbitals, i.e. the intuitive “squared LCAO weight” made rigorous for a non-orthogonal basis (non-negative, summing to exactly 100%). Mulliken’s overlap cross terms instead go negative or overshoot on diffuse empty levels, which is why --projection mulliken is offered but not the default; the two measures agree in kind on occupied levels and disagree strongly on empty ones, where atomic attribution is convention-dependent in the first place.

Symmetry does the orbital selection: a crystal state of irrep Γ only picks up the Γ-adapted combination of each shell (symmetry-forbidden shells project to ~0 and are culled at 0.1%), so every entry carries the crystal irrep — Sc 3d R5+ 76.0%. Projections onto the fragment eigenstates are deliberately not shown as percentages, although they still position the correlation lines and the alignment anchors: a fragment level such as “Sc 3d R5+” itself mixes 3d and 4d AO character, so its weights disagree with the AO populations, and the two lists side by side would read as a contradiction.

The wave-function sketches follow the same measure. Lobe signs and orientations come from the r0 = 2 bohr radial amplitudes, so orthogonalization tails keep their inverted phases, while lobe sizes follow the per-(atom, l) Löwdin population — a diffuse gth-dzvp Sc 4p has 5.4x the amplitude of an F 2p at r0, which would otherwise draw an ~8%-population Sc admixture at 83% of the largest F lobe. The calibration factor of each (atom, l) channel comes from the multiplet-summed populations, so symmetry-equivalent atoms draw exactly equal lobes and the sigma/pi contrast between degenerate partners survives.

Why a valence level can look antibonding (--valence-only)

A valence level can faithfully look “antibonding” around an atom that carries a semicore shell. In ScF3 the F 2p σ R1+ level carries Sc 3s (semicore) c = +0.26 against Sc 4s c = +0.05, so the true wave function — verified against a real-space pbc_eval_gto evaluation — has a radial node 0.7 Å from Sc and is negative near the atom. That is the orthogonality tail against the −50 eV Sc 3s band, not a drawing error.

The scale separation is the ⟨r⟩ one: in this basis ⟨r⟩(Sc 3s) = 0.77 Å against ⟨r⟩(Sc 4s) = 1.86 Å and a 2.03 Å bond — the compact semicore shell cannot bond, and its admixture is on-site orthogonality. --valence-only applies exactly this reasoning to the drawing: shells whose occupied fragment bands lie more than 12 eV below the crystal VBM on the aligned scale (Sc 3s, Sc 3p, F 2s here; detected automatically and printed) are dropped from the drawn wave functions, so the σ level shows its bonding Sc 4s component — Sc s flips from −0.48 to +0.47 against the same F p lobes — while levels a semicore shell dominates (the semicore bands themselves) keep it.

The same asymmetry appears in the COOP colors: semicore orthogonality tails (shells whose own bands lie ≥ 10 eV below their fragment column’s HOMO) are excluded from the overlap population unless the level is that semicore band. With the Sc 3s tail included, the bonding σ R1+ of ScF3 would read P = −0.015 instead of +0.112. Antibonding |P| is systematically larger than bonding |P|, the usual non-orthogonal COOP asymmetry.

What the crystal-orbital regression test asserts

Testsuite section 3 exercises the crystal --diagram --pyscf engine with a deliberately small run, skipped when pyscf is not installed: --diagram --pyscf --onsite on NaCl at X with --kmesh 1 1 1 --ke-cutoff 80, which needs a single crystal SCF and about a minute of wall clock. It asserts the report (one SCF and no fragment SCF, the D+ S D = S verification, the induced representations of every shell, the Cl-3p valence band bonding with Na 3p, the written coupling table) and two invariants of the embedded wave-function sketches, which manual inspection alone used to cover: the drawn σ lobe phases along the Na–Cl bonds must reproduce the bonding/antibonding letter that the engine derives independently from the COOP overlap population, and the Na/Cl ratio of the p-channel lobe sizes must equal the square root of the ratio of the level’s Löwdin populations — the per-(atom, l) lobe-size calibration documented with the command in section 3 of crystod, whose collapse the earlier extended-Hückel tests could not see. A --conventional rerun against the written --chk then covers the conventional-cell display path with no second SCF.

References for the extended-Hückel construction

M. Wolfsberg and L. Helmholz, J. Chem. Phys. 20, 837 (1952); R. Hoffmann, J. Chem. Phys. 39, 1397 (1963).