24The lattice-model route

A deliberate simplification, not automatically first-principles chemistry

Our main nitrogen workflow began with nuclei, a basis set, orbitals, integrals, and a selected correlated model. The repository also supports a complementary route: choose a graph and a small set of interaction parameters, then construct a model Hamiltonian directly.

A model Hamiltonian retains selected physics while deliberately omitting other details. A Hubbard model might retain hopping and local electron repulsion; an Ising model might retain a two-level degree of freedom per site and selected spin couplings. Their parameters can come from experiment, a separate electronic-structure calculation, or an illustrative benchmark. A graph shaped like a molecule is not, by itself, an accurate molecular Hamiltonian.

This branch can bypass AO-basis selection, SCF, and molecular integral evaluation. It does not bypass the need to justify the physics or choose the correct particle/spin sector. We will keep the distinction from neutral singlet N₂ at 1.85 Å explicit: the examples below are alternate models, not cheaper computations of exactly the same nitrogen problem unless a defensible model derivation establishes that relationship.

All construction patterns are unexecuted. Source tests and notebook settings are referenced as records of the repository, not new scientific or hardware results.

Sites, bonds, and boundaries

A site is a labeled degree of freedom. It may represent a localized orbital, a spin, or another two-level subsystem. An edge, or bond, specifies a connection between two sites; its weight can modulate an interaction. A graph stores those sites and connections. Its graph dimension/topology does not automatically specify atomic coordinates, orbital shapes, or interatomic distances.

An open chain has endpoints. A periodic chain adds a wrap bond to form a ring. Periodic boundaries can reduce edge effects in a finite model, but they do not make it an infinite material. Finite-size effects are differences caused by using a limited number of sites; results may change with size and boundary choice. Periodic graph links also do not implement periodic electronic-structure integrals, reciprocal-space sampling, or an Ewald Coulomb sum.

LatticeGraph stores a sparse weighted adjacency matrix, a symmetry flag, site count, and optional edge coloring. It does not store a coordinate array from which the PPP potential routines silently derive distances. Python exposes properties such as num_sites, num_edges, num_nonzeros, and is_symmetric, and methods weight(i,j), are_connected(i,j), adjacency_matrix(), and sparse_adjacency_matrix() through the bindings.

A symmetric adjacency matrix describes reciprocal links. For an ordinary nonzero undirected edge, both directions are stored, so the sparse entry count is twice the unique edge count. num_edges counts upper-triangle entries; it is not a meaningful general directed-edge count. Explicit stored zeros and self-links should not be used to infer physical bonds from these summary counts alone.

Custom graphs can be constructed from an edge-weight dictionary or a dense/sparse square matrix. Dictionary edges are stored as given; they are not automatically made reciprocal. With num_sites=0, the dictionary constructor infers size from the largest index, so an explicit size is needed to retain extra isolated sites. It checks edge indices against that size; matrix constructors check squareness. These are structural checks, not complete validation of physical parameters or finite values.

LatticeGraph.make_bidirectional(graph) computes $A+A^T$. It is useful for a graph supplied in one direction only. Calling it on an already symmetric graph doubles weights rather than merely confirming them. This concrete implementation rule matters because model builders multiply edge weights into their couplings.

The shipped topology constructors

The factory names are chain, square, triangular, honeycomb, and kagome, plus custom matrix/dictionary construction. Their default edge weight is t=1.0, boundaries are open unless requested, and dfs_ordering=False by default. These are data factories, not calls to the chemistry algorithm registry.

A chain has n sites labeled 0 through n−1. periodic=True adds a wrap bond only when n>2; smaller chains are not turned into doubled two-site rings. A square lattice has nx×ny sites with index y*nx+x, connecting right/up neighbors and requested wraps. A triangular lattice adds the upper-right diagonal links, including single-axis and corner wrap cases as appropriate.

Honeycomb uses two sites per unit cell, indexed 2*(y*nx+x) and the following integer. Its intra-cell and inter-cell bonds make the familiar three-neighbor bulk connectivity without storing a drawing of hexagons. Kagome uses three sites per cell, indexed 3*(y*nx+x)+s, and corner-sharing triangle connectivity. These counts are essential: nx×ny honeycomb cells do not mean nx×ny sites.

Zero sizes are rejected. The two-dimensional constructors require an extent of at least two when that direction is periodic. At the smallest periodic sizes, topology generation can emit the same site pair more than once; Eigen triplet assembly adds duplicate weights. For example, a periodic length-two square-lattice direction has both a forward and a wrap contribution on the same pair. Inspect the resulting weights rather than assume every unique edge still has weight t.

Chain, square, and honeycomb factories supply topology-aware edge coloring. Triangular and kagome use randomized greedy coloring with 32 trials and coloring_seed=0 by default. An edge coloring assigns colors so that edges of the same color share no vertex; same-color two-site gates can therefore be scheduled on disjoint sites. This differs from the term-commutation graph coloring of Lesson 20. Open/even rings often use alternating bond colors; odd rings need an additional color. The graph tests verify explicit adjacency, boundary links, and color validity on their fixtures.

For chain, square, and triangular graphs, dfs_ordering=True searches for a path visiting each site exactly once and relabels accordingly. Here Hamiltonian path is a graph-theory name, not an energy operator. The search is depth-first with backtracking and can be expensive; failure raises. Honeycomb and kagome currently ignore the exposed DFS-ordering flag. Do not assume that option changes their ordering.

LatticeGraph.permute(graph,path) returns a permuted adjacency and remapped coloring. Its current implementation assumes a valid complete permutation rather than exhaustively validating it. External site parameters, distances, state labels, and fermionic mode conventions must be permuted consistently; graph reindexing alone does not transform an already mapped fermionic state’s signs.

Two different output routes

The Hückel, Hubbard, and PPP functions return a molecular-style Hamiltonian backed by SparseHamiltonianContainer. Their implementations are native C++ templates in model_hamiltonians.hpp, re-exported to Python by bindings. They supply one-/two-body data that the existing fermion-to-qubit mapper can consume.

The Ising and Heisenberg functions instead directly return QubitOperator. One qubit represents one model spin site, so no additional Jordan–Wigner transformation is required to interpret that spin Hamiltonian. A four-site Hubbard model ordinarily maps to eight spin-orbital qubits; a four-site spin-half Heisenberg model has four qubits. Confusing these two routes changes the Hilbert space.

ModelOrbitals supplies abstract mode counts and active-space indexing without an AO basis, coefficients, or orbital energies. Its implementation defaults to a full active space with an axis-free single channel treated as restricted. Explicit symmetry-blocked index sets can describe other cases. The sparse model container creates this abstract orbital object; its real integral storage is restricted.

Some nominal AO/MO dimension accessors return the model size for compatibility, but that does not imply actual atomic basis functions exist. Accessors needing those data, such as get_basis_set, fail. Population analysis, AO visualization, and molecular geometry workflows cannot be inferred to work merely because the result has an Orbitals-derived object. The ModelOrbitals tests distinguish structural counts from real orbital payloads.

The graph and Hamiltonian builder do not choose electron filling. For fermionic models the populations are supplied to a solver, trial configuration, or sector-selection procedure. The standard mapping still uses both spin modes per site, even though ModelOrbitals’ generic index vocabulary is not itself a declaration of a many-electron spin state.

Hückel and tight binding: moving occupation between sites

For real symmetric hopping parameters, the shipped Hückel form is

Here $n_i=n_{i\alpha}+n_{i\beta}$, $\epsilon_i$ is a site energy, $t_{ij}$ a hopping strength, and $w_{ij}$ the graph weight. “Hopping” moves occupation between modes; it does not require imagining an electron as a classical ball tracing one bond at a time. The model is noninteracting: there is no explicit electron–electron repulsion term in this constructor.

The Python call is create_huckel_hamiltonian(lattice, epsilon=..., t=...). Both physical parameters are required; there are no hidden fitted defaults. Epsilon may be a scalar or a length-n vector; t may be a scalar or an n×n matrix. The code sets diagonal one-body entries to epsilon and off-diagonal linked entries to minus hopping times graph weight. Diagonal graph entries are not used as a substitute for epsilon.

For two linked sites, epsilon=0, graph weight 1, and t=1 give the illustrative one-body matrix

The symmetric orbital has energy −1 and the antisymmetric orbital +1. With two noninteracting opposite-spin electrons occupying the lower orbital, the model energy is −2 before any separately introduced constants. These are hand-worked model numbers, not nitrogen predictions or executed solver output.

A symmetric graph alone does not force a supplied hopping matrix to be symmetric: the implementation uses t(i,j) and t(j,i) separately and checks dimensions, not a universal Hermiticity condition on all parameters. Supply consistent finite real symmetric hopping for this Hermitian model. Graph weight t and model hopping t are separate multipliers; setting both to an energy-valued hopping can accidentally square its intended magnitude. A graph weight of one is a clear baseline.

Hubbard: hopping competes with double occupation

create_hubbard_hamiltonian(lattice, epsilon=..., t=..., U=...) adds

U is required and can be scalar or length n. For positive U, two opposite-spin electrons on the same site pay a repulsion cost. Equal-spin double occupation of one spin orbital is already forbidden by fermionic statistics; this interaction concerns the two different spin modes at a site. Negative U can describe a different attractive model, not ordinary repulsion inferred automatically from a molecule.

The native builder stores $g_{iiii}=U_i$ as raw chemist-order two-body integrals. It does not store U/2: the spin sums and one-half in the second-quantized Hamiltonian already give the correct on-site interaction. The default model scalar is zero. Existing Python tests check hopping signs, on-site integrals, absent intersite integrals, and the core-energy value.

Half filling in this spinful one-orbital-per-site model means n electrons on n sites, since full occupation is 2n. At large positive U compared with hopping, double occupation is suppressed, but quantum mixing through virtual double occupation still affects the low-energy spin state. At small U, delocalization plays a larger role. This is why the ratio U/t is informative; it is not a theorem that quantum computation is advantageous in every large-U case.

A two-site half-filled singlet sector can be expressed in a singly occupied singlet and a symmetric double-occupation combination. With suitable basis phases its illustrative matrix is

whose lower eigenvalue is $(U-\sqrt{U^2+16t^2})/2$ for zero site energy. This explicitly shows that the state mixes configurations and that suppressing double occupation does not mean setting hopping to zero. In the large-U limit the low-energy correction scales like $-4t^2/U$. The repository constructs the Hubbard operator; this analytic reduction is explanatory physics, not a separately invoked effective-spin solver.

A four-site half-filled model with two alpha and two beta electrons has $\binom42^2=36$ determinants in that fixed-population sector, compared with 256 states in its full eight-qubit Fock register. Diagonalizing the full register without a sector restriction can return an energy for another electron count. The earlier classical qubit-solver limitations remain relevant.

PPP and extended density interactions

The Pariser–Parr–Pople model, abbreviated PPP, extends Hubbard with intersite charge interactions:

For symmetric V this is one contribution per unordered pair. The z values are reference/background charges in this model, often one for a one-electron-per-site reference. They are not automatically the physical nuclear charges, an ECP, or frozen-core data from an ab-initio molecule.

create_ppp_hamiltonian requires epsilon, t, U, V, and z. Per-site epsilon/U/z accept scalars or length-n vectors. Per-pair t/V accept scalars or n×n matrices. The hopping still follows graph edges, but the V loop considers all off-diagonal pairs supplied in the matrix, not only graph neighbors and not multiplied by the graph’s hopping weight. A sparse graph does not force a short-range Coulomb interaction.

Expanding the density term explains its storage. The builder adds raw $(ii|jj)=V_{ij}$, modifies the site-energy diagonals for the reference charges, and returns the scalar

as the Hamiltonian core energy. For symmetric V, each diagonal receives $-\sum_{j\ne i}V_{ij}z_j$. This is more than adding a two-body term: omitting the linear and constant parts would describe another reference-charge convention.

In a two-site toy with z=1 and V=0.2, the additional interaction is $0.2(n_0-1)(n_1-1)$. It contributes a +0.2 scalar, −0.2 to each site’s one-body energy, and the intersite density product. When both sites have one electron, those three pieces sum to zero. The constant does not vanish just because the mapped active operator is being sent to a quantum algorithm.

The C++ PPP test makes the all-pairs convention especially clear: a four-site chain with scalar V=1 and z=1, and other parameters zero, has diagonal correction −3 per site and offset 6. Those are recorded test expectations, not a run here. The chain’s three graph bonds did not limit V to three pairs.

Use a symmetric physically justified V. Parameter shape validation does not supply a complete physical symmetry/positivity proof, and the sparse mapper’s symmetry handling is not a repair procedure for arbitrary inconsistent tensors. Ordinary qubit mapping excludes the stored core scalar, including this PPP background offset; add it back once when reporting the model’s full energy.

Screening functions need supplied distances

The helper pairwise_potential(lattice,U,R,func,nearest_neighbor_only=False) evaluates unique pairs, uses $U_{ij}=\sqrt{U_iU_j}$, reads the supplied R(i,j), calls the user formula, and stores a symmetric result. Diagonal entries remain zero. Scalar R assigns the same distance to every evaluated pair; it does not calculate graph shortest-path or Euclidean distances.

The shipped Ohno and Mataga–Nishimoto formulas are

They use Hartree for U, Bohr for R, and dimensionless epsilon_r=1.0 by default. For positive physical parameters they approach the on-site scale at small separation and screened inverse distance at large separation. Setting nearest_neighbor_only=True limits evaluation to connected pairs; the default is all pairs. These behaviors come from the native potential functions, not an automatic dielectric calculation.

The code validates parameter shapes but does not determine a material’s permittivity, repair unphysical distances, or guarantee a real geometric mean for arbitrary signed U inputs. A phenomenological screening law can be a useful approximation and still require calibration against the question being asked. Periodic graph boundaries do not add periodic electrostatic images to these formulas.

Ising and Heisenberg: direct spin Hamiltonians

The Python spin-model utilities build Pauli operators directly. For create_ising_hamiltonian(graph,j,h=0.0), the convention is

J is required; h is a transverse field with default zero. There is no hidden leading minus sign. Positive J in this convention favors opposite Z values on a single bond; a textbook ferromagnetic model written as $-JZ_iZ_j$ with positive J requires a negative supplied coupling. The direct Ising helper does not expose an independent longitudinal-Z-field parameter.

create_heisenberg_hamiltonian requires jx, jy, and jz, with hx, hy, hz defaulting to zero. It builds

This includes anisotropic exchange and all three local field directions. Equal exchange components give the isotropic case. It does not include arbitrary cross-axis or antisymmetric exchange terms merely because the general Pauli data class could represent them. An Ising model with a longitudinal field can be expressed through the appropriate Heisenberg parameter choices rather than a nonexistent Ising keyword.

The formulas use Pauli matrices. If a physics convention writes $J\mathbf S_i\cdot\mathbf S_j$ for spin one-half with $\mathbf S=\boldsymbol\sigma/2$, pass exchange coefficients J/4 in this builder. A field coefficient multiplying S rather than sigma likewise differs by two. These factors are essential when transferring parameters from a paper or experiment.

Pair couplings accept a scalar or n×n matrix; fields accept a scalar or length-n vector. Each undirected edge is counted once. Spin-model pair construction reads the upper-triangle coupling for that edge; it does not average inconsistent upper/lower inputs. Graph weights multiply bond couplings but not on-site fields. For a two-site toy with J=1 and zero transverse field, ZZ gives energy +1 on aligned bits and −1 on opposite bits. These are spin configurations, not electron occupation sectors of a four-mode Hubbard problem.

Geometry-aware scheduling is metadata, not new physics

include_term_groups=True is the default for Ising/Heisenberg. When graph edge coloring is present, the constructor creates a LayeredPartition: field-direction groups have disjoint one-site terms, and each exchange-axis group is divided into edge-color layers. Coefficients and labels still describe the same operator. The alternate ungrouped path uses symbolic Pauli arithmetic and converts canonical strings into the opposite QubitOperator label direction explicitly.

A custom graph has no automatically stored topology coloring. The constructor then falls back to ungrouped construction rather than discovering every possible optimum. The free coloring functions and optional external NetworkX helpers offer other tools, but graph coloring and state-dependent measurement optimization are different tasks.

The layered partition can help an evolution builder preserve a useful schedule. The current expectation estimator does not jointly measure that hierarchy; it falls back to individual terms unless supplied with a compatible flat QWC partition. A geometric gate-layer advantage is not automatically a shot-count advantage. The model tests verify terms, nonuniform couplings, and graph-weight multiplication, while the partition tests discussed earlier cover downstream grouping behavior.

The cyclobutadiene notebook: chemistry-informed parameterization

The source extended_hubbard.ipynb motivates a four-site model from a square cyclobutadiene pi system. It first uses molecular orbitals, optional PySCF AVAS selection targeting carbon 2pz functions, and native Pipek–Mezey localization to obtain site-like orbitals. Thus this particular model workflow does not bypass all molecular preprocessing; it uses that preprocessing to motivate a simpler parameterized operator.

The notebook reconstructs the Hamiltonian from localized orbitals, takes the largest absolute off-diagonal one-body entry in each row as a nearest-neighbor hopping indicator, averages those values into t, and averages diagonal two-electron entries into U. It then assumes a periodic four-site chain, averages designated C–C distances from the structure in Bohr, and generates nearest-neighbor-only V with the Ohno formula. Epsilon is set to zero and z to one.

Those choices compress a much richer integral tensor into a few parameters. They discard heterogeneity, hopping signs/phases and many interaction channels; extracting magnitudes is not a general proof that the resulting ring is unitarily equivalent to the ab-initio active Hamiltonian. Localized orbital ordering is not automatically atom ordering. The measured/averaged bond-distance assumptions also rely on the input atom order. Parameter variation should assess those choices rather than treat U/t as an automatic industrial-validity or quantum-advantage metric.

Setting a uniform epsilon to zero changes energy by a constant only within a fixed electron-number sector, where $\epsilon\sum_i n_i=\epsilon N$. It is not a global identity shift across the full Fock space. Comparing absolute model energies or allowing another filling requires remembering that omission.

The notebook then solves the small model using multi_configuration_calculator/macis_cas with explicitly supplied active alpha/beta counts, prepares selected trial determinants, maps to qubits, and studies IQPE. A localization output is a determinant carrier, not automatically a transformed correlated CI state; the fresh model solve is essential. The optional PySCF selection and the native localizer’s limits from earlier lessons remain applicable.

Its descriptions of molecular instability or reaction pathways motivate the exercise, but a parameterized square-ring calculation does not alone establish a quantitative automerization barrier, room-temperature lifetime, or useful chemical process. Those require matched geometries, energetics, and additional physical corrections. No notebook cells or output-generating operations were run here.

A larger Ising example is a circuit study, not a statevector simulation

The source estimation_ising_2d.ipynb chooses an open 10×10 square graph: 100 sites and 180 unique nearest-neighbor bonds. It supplies ZZ coupling 1 and X field 0.5, giving 100 problem qubits. The declared time-independent drive is one; its evolution builder uses one outer interval, fourth-order Trotter, and two inner divisions, starting from all-zero spin state.

The arithmetic operator inventory is 180 bond terms and 100 field terms before schedule expansion. A generic 100-qubit statevector would have $2^{100}$ complex amplitudes, but constructing/counting a circuit does not require storing that full vector. The notebook compares geometry-grouped and ungrouped circuit resource studies. It is not evidence that this session simulated 100 interacting spins or demonstrated a device speedup.

A printed physical-resource snapshot depends on the selected logical workload, synthesis assumptions, QEC model, and hardware parameters. We will analyze those inputs in Lessons 27–29 rather than treating the graph size or notebook output as a universal feasibility claim.

A small route you can carry into the rest of the course

The following unexecuted pattern follows the native utility and mapper contracts:

from qdk_chemistry.algorithms import create
from qdk_chemistry.data import LatticeGraph, MajoranaMapping
from qdk_chemistry.utils.model_hamiltonians import create_hubbard_hamiltonian

graph = LatticeGraph.chain(2)  # unit graph weight, open boundary
model = create_hubbard_hamiltonian(graph, epsilon=0.0, t=1.0, U=4.0)
mapping = MajoranaMapping.jordan_wigner(num_modes=2 * graph.num_sites)
qubit_model = create("qubit_mapper", "qdk").run(model, mapping)
# A small fixed-population reference could use:
# energy, state = create("multi_configuration_calculator", "macis_cas").run(model, 1, 1)

The model parameters are illustrative consistent energy units, not fitted nitrogen values. Filling is separate from graph construction. After mapping, one can use the same state preparation, time-evolution, measurement, and QPE layers already studied, but only under their compatibility and statistical constraints. The ordinary qubit mapper still excludes model core/background constants; small exact references must use the same sector and offset.

Long-range density interactions can add quadratic numbers of pair terms even when hopping connectivity is sparse. Increasing site count grows fermionic or spin Hilbert spaces exponentially, while the input Hamiltonian can remain compact. Model sparsity, low rank, graph locality, and easy classical solvability are different properties. Special structures may make a benchmark classically tractable, while frustration or correlation may make others difficult; neither fact proves the best quantum implementation wins in practice.

The lattice route is therefore valuable both for interpretable physics and for controlled algorithm benchmarks. It lets us vary connectivity, filling, signs, interaction strengths, and boundaries while keeping the computational construction explicit. Its scientific credibility comes from validating those choices against the intended system—not from borrowing the class name Hamiltonian from the ab-initio workflow. The next lesson examines what it actually means to execute the resulting circuits and connect them to external frameworks.