08Choosing more useful orbitals

The same orbital space can be described in different ways

Our original N₂ reference has 28 spatial orbitals in cc-pVDZ. SCF naturally supplies canonical orbitals, which diagonalize the effective mean-field operator. These are convenient for orbital energies and perturbative formulas, but they need not be the most useful orbitals for describing bonds, comparing geometries, selecting a reduced model, or constructing compact approximate wavefunctions.

This lesson changes the orbital description, not the geometry. The nitrogen nuclei remain 1.85 Å apart unless an explicitly separate example says otherwise. We will distinguish three goals: spatial localization, diagonalization of one-electron density, and optimization of more specialized representation objectives. A method called a localizer in the registry can serve any of these goals; not all of its outputs are spatially localized.

No localization, density calculation, cube evaluation, or test was executed during preparation. Code excerpts describe supported calls and their source-level limitations. Illustrative matrices are mathematical examples, not nitrogen outputs.

What an orbital rotation does—and what it does not preserve

Let the columns of $C$ contain molecular-orbital coefficients in the AO basis, and let $S$ be the AO overlap. Starting from $C^{\mathsf T}SC=I$, a real orthogonal rotation $U$ gives

The rotated orbitals span the same selected space. They are not new basis functions fitted to the atom, and the AO basis size has not decreased. Consequently, a rotation alone does not reduce the number of electron spin modes or qubits in an unreduced representation.

If $U$ mixes only occupied orbitals of a determinant, that determinant changes only by the determinant of $U$, an unobservable overall phase for a unitary transformation. Rotations among unoccupied orbitals leave the occupied determinant unchanged as well. In UHF this applies separately within the occupied alpha and beta spaces. Mixing occupied and unoccupied orbitals generally changes the determinant and its energy; that was the purpose of the stability perturbations in Lesson 6.

A fully represented many-electron state can be expressed in another orbital basis without changing the physics, but its determinant amplitudes must be transformed too. Likewise, an exact complete-active-space calculation is invariant under rotations inside the same complete active space when its Hamiltonian and wavefunction are treated consistently. An active space is the subset of orbitals retained for explicit many-electron correlation; inactive orbitals are fixed occupied orbitals, while excluded orbitals lie outside that correlated model. Mixing across those boundaries changes a truncated model unless its definition is changed consistently.

For a one-particle density matrix $\gamma$ in the MO basis, consistent re-expression is

The AO density, and therefore a consistently defined Mulliken population, is unchanged. Keeping old determinant coefficients after changing occupied–virtual orbitals is not this transformation. Neither is retaining old active-space labels while mixing their orbitals with an excluded space.

The QDK localizer return contract is especially important

The public Python class is OrbitalLocalizer, backed by the C++ Localizer interface. The call is

transformed_reference = localizer.run(wavefunction, loc_indices_a, loc_indices_b)
transformed_orbitals = transformed_reference.get_orbitals()

The result is a Wavefunction, not an Orbitals object despite some older docstrings. The two index lists select columns to transform. Use sorted, unique, valid zero-based indices. For a restricted common spatial-orbital set the lists must agree. An empty list means no orbitals in that channel, not “use all.”

Crucially, current localizers return a single Aufbau determinant carrier with updated orbitals. They do not in general transform and preserve an input CI expansion or coupled-cluster amplitudes. Some methods attach derived, transformed 1-RDM data, but this does not turn the carrier’s single determinant into the original correlated state. Even an empty selection can return such a new single-reference carrier with unchanged orbital coefficients. The factory helpers and implementation make that behavior explicit.

Thus a correlated input should be retained. Use the returned orbitals to build a new Hamiltonian and solve the chosen correlated problem in that representation before interpreting new CI coefficients, energies, or entropies. The fact that a transformed 1-RDM is attached is not proof that every higher-order property of the original state survived.

At this revision the native registry provides qdk_pipek_mezey, qdk_vvhv, qdk_natural_orbitals, qdk_gauge_fixing, qdk_active_space_qio, and the deprecated qdk_mp2_natural_orbitals. The optional PySCF multi-method implementation is named pyscf_multi, not pyscf. Use available("orbital_localizer") in an ordinary installed environment to discover which are registered. This course has inspected the registrations rather than importing the package.

Canonical orbitals solve a different optimization problem

Canonical SCF orbitals diagonalize the effective Fock or Kohn–Sham operator. Localization instead chooses a different basis inside an allowed subspace by optimizing a localization objective. A lower localization objective is not a lower HF total energy; for occupied-only rotations the determinant energy is unchanged to numerical accuracy.

A canonical orbital may extend over both nitrogen atoms and have definite molecular symmetry. That is not a flaw. Delocalization can be the cleanest description of the mean-field operator. Localized combinations can instead make bond and atom-centered character easier to see. Degenerate or nearly degenerate orbital energies also make individual canonical columns nonunique or sensitive to small changes in geometry.

After a nontrivial rotation, the original list of orbital energies does not generally diagonalize the transformed mean-field matrix. Native Pipek–Mezey, natural-orbital, and VVHV outputs invalidate orbital energies. The current PySCF localizer adapter copies input energy arrays into its output, so their presence must not be interpreted as proof that the localized columns are canonical. Algorithms requiring canonical orbital energies need a scientifically consistent canonical reference, not merely has_energies()==True.

Pipek–Mezey concentrates an orbital’s atomic populations

Pipek–Mezey (PM) localization favors orbitals whose population is concentrated on a smaller number of atoms. For one normalized orbital $p$, define its Mulliken population on atom $A$ as

This is the population of one orbital with unit normalization, not a net atomic partial charge and not a factor-of-two closed-shell electron count. Native PM maximizes

For illustration, an orbital split equally between two atoms contributes $0.5^2+0.5^2=0.5$; an orbital wholly assigned to one atom contributes one. Maximizing the sum therefore favors concentration under the chosen population rule. Bond orbitals can still extend across two atoms; “localized” does not mean each orbital becomes a delta function on a nucleus.

The native implementation uses Jacobi, or pairwise plane, rotations. For two orbitals,

It computes atomic populations and cross-populations, derives candidate rotation angles and objective increases for orbital pairs, applies advantageous rotations, and repeats sweeps. The overlap-times-coefficient matrix is rotated alongside the coefficients, avoiding a full recomputation from scratch for each pair.

The internal stopping controls have defaults tolerance=1e-6, max_iterations=10000, and small_rotation_tolerance=1e-12; the schema requires at least two iterations. tolerance is used for objective-change and small proposed angle/increase tests, not an electronic energy accuracy. small_rotation_tolerance zeros small intermediate pair coefficients before choosing angles. The internal routine logs whether it converged, but returns its last coefficient matrix rather than throwing solely because the sweep limit was reached.

There is a concrete current limitation: the public qdk_pipek_mezey wrapper constructs a fresh internal default settings object rather than forwarding its own configured settings to the iterative routine. Those exposed settings can be stored, but changing them is not honored by this wrapper’s calculation at this revision. We therefore do not present a changed tolerance as a verified control of native PM convergence. The VVHV path below does forward its settings to its internal PM object. This is a source-level behavior to account for, not a request to modify the repository.

Native PM requires the molecular basis, its atom mapping, and AO overlap. It localizes the selected columns as one group. It does not automatically split a mixed input list into occupied, inactive, active, or virtual subsets. If preserving the HF determinant is the goal, select occupied and virtual groups separately. If preserving a later frozen-core partition is also important, separate core and valence occupied groups too.

For our original restricted HF reference, localizing only its seven occupied orbitals is a valid determinant-preserving illustration:

from qdk_chemistry.algorithms import create

n_alpha, n_beta = reference.get_total_num_electrons()
occupied = list(range(n_alpha))
pm = create("orbital_localizer", "qdk_pipek_mezey")
occupied_localized = pm.run(reference, occupied, occupied)

This assumes the original RHF reference, not an unrestricted result from the stability workflow. It changes occupied orbital shapes, not the electron count or molecular geometry. The native localization tests explicitly check metric values, overlap orthonormality, selected-space unitary transformations, and separate occupied/virtual calls.

Boys, Edmiston–Ruedenberg, and Cholesky offer different meanings of local

The optional implementation is configured, for example, as

boys = create("orbital_localizer", "pyscf_multi", method="foster-boys")
boys_reference = boys.run(reference, occupied, occupied)

The adapter dispatches to PySCF’s lo.PM, lo.Boys, lo.ER, or lo.cholesky_mos. QDK supplies the converted basis and selected coefficient columns; the external library supplies those localization implementations.

Foster–Boys minimizes total orbital spatial spread:

The first term measures the squared distance distribution and the second removes the squared position of its center. For rotations within a fixed subspace, the sum of the first terms is invariant, so one can equivalently maximize the sum of squared orbital centers. Dipole/moment integrals connect this spatial objective to the coefficient matrices. This can produce compact bond or lone-pair descriptions without relying on Mulliken atom assignments.

Edmiston–Ruedenberg (ER) instead maximizes the sum of each orbital’s Coulomb self-repulsion:

Concentrating an orbital’s density tends to increase that self-repulsion, making it a localization objective. We are maximizing this auxiliary quantity while not changing the physical occupied determinant under an occupied-only rotation. This is not a proposal to maximize the molecule’s total electron-repulsion energy. ER requires two-electron integral work and is generally more expensive than population- or moment-based localization.

Cholesky localization, selected by method="cholesky", is a noniterative alternative based on a rank-revealing factorization of the selected orbital subspace’s coefficient-density/projector representation. Pivoting identifies independent AO-associated directions, and factorization constructs a localized spanning set. Conceptually one factors a matrix built from $C_{\mathrm{sel}}C_{\mathrm{sel}}^{\mathsf T}$; in the complete-rank case the new factors span the same selected subspace. It is cheaper than optimizing every orbital-pair objective repeatedly, but does not guarantee the optimum of Boys, PM, or ER. This is not Cholesky factorization of the four-index electron-repulsion tensor, which is a different subject in Lesson 12.

For PySCF PM only, population_method defaults to "mulliken" and accepts "lowdin", "meta_lowdin", and "becke" as well. Löwdin analysis assigns populations after symmetric AO orthogonalization; meta-Löwdin is an external modified orthogonalization/atomic-partition scheme; Becke uses smooth atom-centered real-space weights rather than assigning overlap terms solely by AO labels. These change the PM localization criterion. They do not add those methods to the standalone population_analyzer, which exposed only Mulliken in Lesson 7.

The adapter defines occupation_threshold=1e-10, but its current localization body does not use it to partition the chosen columns. Despite broad docstrings about automatic separation, callers must choose their groups. Its schema also does not forward arbitrary PySCF localization iteration tolerances; external defaults govern the kernel calls shown here.

A small but important input limitation is visible in the current code: the helper’s single-index path returns the full input matrix before the caller copies its first column back to the selected slot. A nonzero singleton selection can therefore copy column zero into that slot rather than act as a no-op. Do not use singleton calls to this adapter as an assumed harmless operation; skip localization of one-column groups at the calling level. This is particularly relevant when separately treating one singly occupied ROHF orbital. Multi-column occupied selections such as the nitrogen example follow the actual localization kernels.

Natural orbitals diagonalize density, not spatial spread

The 1-RDM introduced in Lesson 7 can be written in many orbital bases. Natural orbitals diagonalize it:

$n$ is diagonal, and its entries are natural-orbital occupation numbers, often abbreviated NOONs. For a physical electron 1-RDM in a common spatial-orbital basis, spin-traced occupations lie between zero and two; spin-resolved occupations lie between zero and one. Their sum is the corresponding electron count. Natural orbitals need not be spatially localized.

An illustrative matrix

has occupations 1.3 and 0.7 in the equal-weight plus/minus orbital combinations. Off-diagonal coherence in the original basis becomes unequal diagonal populations in the natural basis. This is a density-matrix example, not the two-determinant Hamiltonian from Lesson 2 and not a computed nitrogen density.

For an RHF determinant, the spin-traced 1-RDM has occupations two and zero, so the occupied and virtual natural-orbital subspaces are highly degenerate. Natural-orbital analysis of that determinant alone cannot reveal the correlated occupations it never contained. Fractional occupations from a correlated density can signal that a single paired determinant is inadequate, but the interpretation depends on where the density came from.

In particular, a spin-traced UHF/UKS density can have fractional spatial occupations even though its full spin-orbital state is one determinant. Adding two different spin projectors need not give a spatial matrix with only zero/two eigenvalues. Thus fractional spatial NOONs are not by themselves proof that an exact correlated wavefunction has been obtained. Approximate perturbative density matrices also require physical-consistency checks.

The native qdk_natural_orbitals implementation has no configurable algorithm settings. Its actual requirements are stricter than “any occupation vector”: matching alpha/beta active and inactive spaces, identical selected lists exactly equal to the active-space indices, AO overlap, and real spin-dependent active 1-RDM blocks of the correct dimensions. It adds the spin blocks, symmetrizes small numerical asymmetry, diagonalizes the active density, and sorts occupations in descending order. Only the active coefficient columns are replaced.

For unrestricted input, it reconstructs each active spin density in the AO basis and projects both into the selected alpha spatial-orbital span before diagonalization. The output is restricted. If the beta active span is not contained in that target span, this projection is not a guaranteed lossless re-expression of the whole unrestricted active density. Full common-space examples avoid that problem, but matching orbital indices alone does not prove matching spans.

The output invalidates canonical energies and carries rotated spin-resolved and spin-traced 1-RDM data, together with its single-reference carrier. Preserve the original correlated state and re-solve for new correlated amplitudes when needed. The natural-orbital regression uses a deliberately spin-broken PBE/UKS N₂ calculation at 2.0 Å and verifies descending fractional occupations in the full space. That is a test of this capability, not a result for our 1.85-Å HF calculation.

The older MP2-natural-orbital shortcut is still registered, but deprecated

qdk_mp2_natural_orbitals constructs a selected-space Hamiltonian and asks MACIS for a density-derived orbital transformation from second-order Møller–Plesset perturbation theory, MP2. MP2 estimates correlation through excitations from a canonical reference, weighted by electron interactions and orbital-energy differences. Near small denominators, such as strongly stretched-bond situations, a perturbative approximation can become unreliable; its natural orbitals are not automatically superior.

The legacy implementation requires restricted closed-shell input with orbital energies, matching sorted selections, and both occupied and virtual orbitals in the chosen subset. It constructs integrals, calls macis::mp2_natural_orbitals, rotates the selected columns, clears orbital energies, and attaches a spin-traced occupation-density payload. Electron counts come from the input wavefunction, not missing user settings suggested by an old error message.

Its deprecation guidance is to use qdk_natural_orbitals with a wavefunction that already supplies the MP2 active-space 1-RDM. That separates generation of an approximate correlated density from the generic diagonalization step. A suitable MP2 RDM must be generated or imported separately. The current native dynamical-correlation MP2 calculator returns amplitude storage, not an RDM-bearing state; Lesson 11 explains that boundary.

VVHV separates useful virtual character from excess flexibility

Virtual orbitals are unoccupied in the reference determinant, not physically irrelevant. They provide directions for excitations and orbital response. In a large flexible basis, however, ordinary localization of the entire virtual space can mix chemically interpretable valence-antibonding directions with many high-flexibility functions. The objective may become poorly tied to the chemical distinction we care about.

Valence-virtual/hard-virtual separation (VVHV) uses a minimal basis as a reference vocabulary. Valence virtuals are the part of that compact vocabulary orthogonal to occupied orbitals. Hard virtuals complete the remaining full-basis virtual space. “Hard” is a construction label for extra basis flexibility, not a theorem that these orbitals never matter for chemistry or correlation.

The native implementation is "qdk_vvhv". Its current call contract requires all virtual orbital indices, not all molecular orbitals, in each applicable spin channel. Some header text says all orbitals; the implementation and successful tests pass only the virtual range. It leaves occupied coefficients unchanged. For a restricted reference, the two lists must match; for unrestricted input it constructs each spin’s virtual complement separately. It currently rejects a reduced-rank basis where AO and MO counts differ.

Let $S_{FF}$ be overlap in the full basis and $S_{FM}$ cross-overlap between full and minimal basis functions. The construction first represents the minimal functions in the full AO basis by solving a least-squares problem equivalent, in a well-conditioned full-rank case, to

It orthonormalizes this compact representation in the full-basis metric. With occupied coefficients $C_o$, it projects away occupied components:

After removing the expected near-null directions and orthonormalizing the remainder, these columns form valence virtuals. This explains the linear algebra: subtract components already represented by occupied orbitals, rather than deciding valence character from an arbitrary virtual energy cutoff.

The current inner iterative localizer is Pipek–Mezey, explicitly constructed in the wrapper. It localizes the valence-virtual set. Although an internal class accepts an inner-localizer object, the public settings do not expose a choice of Boys or ER for that role. Occupied localization is a separate call before VVHV, not something the present VVHV wrapper automatically performs.

Hard virtuals are constructed atom by atom and angular-momentum block by angular-momentum block. The implementation builds prototype hard directions, projects atom-centered functions out of the combined occupied-plus-valence-virtual space, aligns the remaining directions with the prototypes, and orthogonalizes them globally. This supplies an AO-linked description of the full complement rather than a freely rotating virtual set with no chemical reference.

minimal_basis defaults to "sto-3g" and accepts "sto-3g*". weighted_orthogonalization=True computes orbital spreads from dipole and quadrupole integrals and weights each hard-virtual column by the inverse of its spread plus 1e-6 before the final orthogonalization. This biases the resulting orthogonal set toward compact directions. The inherited iterative tolerances are forwarded to the inner PM routine on this path. Projection orthogonalizations also use fixed numerical cutoffs and check expected rank separation; some separation problems are warnings, while size/mapping/rank incompatibilities are errors. Neither a completion log nor a successful return proves that a chosen minimal basis defines the best scientific partition.

For nitrogen, the stored STO-3G resource contains two s functions and one p shell per atom: five spatial functions per N, ten for N₂. With seven occupied orbitals, the ideal full-rank count is three valence virtuals and $28-10=18$ hard virtuals, totaling the original 21 virtuals. These are resource-derived counts, not a new localization result. No spin mode has been removed.

A correctly scoped unexecuted call is:

n_mo = reference.get_orbitals().get_num_molecular_orbitals()
virtual = list(range(n_alpha, n_mo))
vvhv = create("orbital_localizer", "qdk_vvhv", minimal_basis="sto-3g")
partitioned_reference = vvhv.run(occupied_localized, virtual, virtual)

The VVHV tests use this occupied-PM-then-virtual-VVHV pattern and test active-space metadata preservation. Preserved metadata means the same indices are stored; if the chosen virtual transformation crosses a later active/excluded boundary, that alone does not preserve the chemical content of a truncated model. Plan the partition and transformation together.

Natural-orbital degeneracy leaves a gauge choice

A gauge choice here means a choice of orbital basis inside a subspace that has no unique individual vectors under the preceding criterion. If several natural occupations are equal, any orthogonal rotation within that block keeps the spin-traced 1-RDM diagonal. The subspace is determined, but its individual orbitals are not.

The native qdk_gauge_fixing method makes that choice more reproducible and then targets a quantum-representation objective. It requires restricted natural orbitals, real nearly diagonal active 1-RDM data, overlap, and identical valid selected lists. A selection may cover a subset of the active space but cannot split an occupation-degenerate block.

The implementation first anchors selected blocks to independent AO-associated directions. It then performs coordinate-descent sweeps over pair rotations within those blocks. Each trial is judged by a mapped coefficient norm

where $h_\ell$ are coefficients of the mapped active-space Hamiltonian in a qubit-operator expansion. We will derive that expansion later; for now, this is a measure of the size of the operator description, not an electronic energy. The current objective sets the separate core-energy offset to zero and is evaluated before qubit tapering.

The native code builds a reference Hamiltonian once after anchoring, transforms its one- and two-electron integrals for candidate rotations using MACIS utilities, maps them with the native Jordan–Wigner machinery, and sums coefficient magnitudes. Each objective evaluation is therefore substantially more expensive than a small matrix norm. It uses angle_samples=32 over $[0,\pi)$, followed by scalar golden-section refinement, and accepts rotations only when the objective improves by more than improvement_tolerance=1e-10 Hartree.

degeneracy_tolerance=1e-6 controls the allowed occupation spread in each block. max_sweeps=3 limits coordinate-descent passes; zero requests AO anchoring only. hamiltonian_constructor is a nested reference, native qdk by default. The settings enforce at least four angle samples, nonnegative sweep count, and finite meaningful tolerances through schema and runtime checks.

The initial AO anchoring need not improve $\lambda$ relative to an arbitrary input orientation; the monotonic acceptance rule applies to the following optimization from that anchor. Nor does reducing this one objective prove a lower complete resource estimate: synthesis precision, state preparation, tapering, repetitions, and the chosen quantum algorithm also matter. The method carries transformed 1-RDM information but still does not carry the original CI expansion. The gauge-fixing tests verify occupation spectra and CASCI energy invariance by rebuilding and re-solving, rather than treating the returned determinant carrier as the correlated answer.

Quantum-information orbitals use local uncertainty as an objective

The remaining native method, qdk_active_space_qio, rotates a fixed restricted active space to minimize the sum of single-orbital entropies. QIO means quantum-information orbitals. It is not an active-space selector and does not implement full-space QIO or a method that mixes across the active boundary.

For a spatial orbital, the local occupations are empty, alpha-only, beta-only, and doubly occupied. Let $n_\alpha,n_\beta$ be its spin populations and $d$ its double-occupation expectation, which requires two-particle information. The four probabilities are

The local entropy is $s=-\sum_k p_k\ln p_k$, with zero-probability contributions taken as zero. The shared entropy helper uses natural logarithms. The objective is $\sum_i s_i$ over active orbitals. A rotation changes which orbital subsystems carry uncertainty and correlation; it need not reduce the overall many-electron correlation or localize orbitals in physical space.

The QIO implementation requires real spin-dependent active 1- and 2-RDMs, restricted orbitals with overlap, and matching selected indices exactly equal to the active space. Restricted open-shell input is supported; unrestricted orbitals are rejected. It tests pair rotations using a coarse and then fine angle scan, accepts entropy-reducing moves, rotates working RDM tensors, and accumulates one overall orbital rotation.

Its actual controls are max_cycles=200, convergence_tolerance=1e-10 for sweep-to-sweep entropy-sum change, coarse_angle_step=0.02 radians, fine_samples=201, and improvement_tolerance=1e-12 for accepting a pair improvement. The coarse step is bounded between 1e-4 and $\pi/2$, fine samples must be at least four, and the tolerance settings are checked for finiteness. Reaching a cycle limit returns the accumulated transformation; it is not a global entropy-minimum certificate.

Two-particle density storage grows as the fourth power of active orbital count, and repeatedly rotating/inspecting it is not free. The method keeps the input physical state’s density information fixed for the representation search; it does not solve a new correlated eigenproblem after every trial angle. Output contains the transformed orbitals and rotated 1-RDM payloads, generally not diagonal, but not the original correlated coefficients or a retained transformed 2-RDM payload. The QIO tests explicitly run CASCI again before comparing final entropy and energy. That is the correct conceptual handoff, and the full RDM machinery will be taught in Lesson 9.

See the orbital without mistaking the picture for an observable

A cube grid samples a scalar field in a rectangular region of space. For an orbital it evaluates $\phi_p(\mathbf r)=\sum_\mu C_{\mu p}\chi_\mu(\mathbf r)$; for a density it evaluates the corresponding quadratic expression. An orbital isosurface is a surface at a chosen positive or negative amplitude. Its colors commonly distinguish sign/phase, not positive versus negative electron charge. Changing an orbital’s overall sign changes the picture’s colors but not its one-orbital density.

The native CubeGrid and CubeGenerator use Bohr for origin, spacing, and margin. CubeGrid.from_basis_set encloses the nuclear bounding box plus a margin; spacing divides each padded extent by the number of intervals. The default grid is 80 points per axis with a 3-Bohr margin. A one-point axis has zero spacing. Axes are Cartesian-aligned, and the backend imposes an INT_MAX total-point limit; that limit is not a practical memory recommendation.

The generator translates the basis once and uses GauXC’s host orbital evaluator, a different role from its DFT quadrature routines. It supports correctly represented spherical and Cartesian Gaussian basis data; that evaluator capability does not remove the earlier limitations of high-level SCF/adapter Cartesian workflows. Coefficient ordering and normalization must match the basis. For ECP models it samples the explicitly represented density; removed core density is not reconstructed.

A complete no-output-file API illustration for one localized orbital is:

from qdk_chemistry.data.symmetry import SymmetryLabel, axes
from qdk_chemistry.utils import CubeGenerator, CubeGrid

local_orbitals = occupied_localized.get_orbitals()
basis = local_orbitals.get_basis_set()
alpha = SymmetryLabel([axes.alpha()])
coefficients = local_orbitals.coefficients().block([alpha, alpha])
grid = CubeGrid.from_basis_set(basis, 40, 40, 40, 3.0)
field = CubeGenerator(basis).orbital(coefficients[:, 6], grid)

The Python binding returns an array shaped (40, 40, 40) here. Column six is a selected localized orbital; after localization it should not automatically be called the highest-energy occupied canonical orbital. No grid evaluation was performed in preparing this lesson, even though this particular lower-level API can return a field without writing a cube file.

CubeGenerator.density uses the supplied AO density exactly as given. It does not add a missing spin factor or sum separate spin blocks. Supplying half the restricted total density produces half the density field with exactly the same array shape. The cube tests compare AO evaluation with PySCF, test that a rank-one density gives the squared orbital, and verify linear scaling of the supplied density matrix.

For file-oriented workflows, generate_orbital_cubes writes one zero-indexed file per requested spatial orbital, with separate alpha/beta files for unrestricted orbitals. The Python convenience helper generate_cubefiles_from_orbitals defaults to a 40-cubed grid and can use native or optional PySCF evaluation. Its output_folder=None mode returns cube text, but still creates and removes intermediate files; it is not a disk-free API. It was not called under this task’s read-only constraints. Full viewer/UI integration is reserved for Lesson 32.

Choose orbitals for a purpose and check the relevant invariants

For the current nitrogen story, occupied PM orbitals can help interpret the mean-field reference; VVHV can organize the remaining virtual flexibility without throwing it away. A correlated 1-RDM can later produce natural occupations for evidence-based active-space decisions. Gauge fixing and QIO optimize still different objectives. No single orbital set is guaranteed best for every subsequent solver or quantum algorithm.

Useful checks depend on the claim. For a purported rotation within a selected space, verify overlap orthonormality and the selected-space transformation. For occupied-only HF localization, compare the occupied AO density and energy rather than individual column values, which may differ by signs and rotations. For an active-space transformation, preserve its boundaries and re-solve the correlated model before comparing CI energies. For visualization, check grid units, spin factors, and basis ordering before trusting attractive lobes.

Along a bond-stretching scan, orbitals can change sign, exchange order, or rotate inside degenerate subspaces. Localized or AO-anchored descriptions can help maintain chemical interpretation, but the current per-geometry localizers do not guarantee a globally smooth matching of orbitals across geometries. VVHV supplies a minimal-basis reference, not an automatic proof that orbital index twelve denotes the same chemical direction everywhere. Natural-orbital degeneracies are another reason to compare subspaces and character rather than index alone.

Finally, locality can reduce the importance of some long-range integrals or make selected CI amplitudes more compact, but Coulomb interactions remain long ranged and improvements are conditional. A different orbital representation can also make an approximation less compact. Claims of fewer terms, easier state preparation, or lower resource estimates must be tested for the chosen downstream algorithm, truncation, and precision—not inferred merely from the word “localized.”

We now have the vocabulary to separate basis functions, molecular-orbital coefficients, occupation/density information, and determinant amplitudes. The next lesson turns that last distinction into an actual correlated electronic solution: configuration interaction, selected determinant methods, and the reduced density matrices that many of these orbital transformations require.