API Reference

rotation2D(θ)

Returns the 2x2 rotation matrix R = R^T parametrized by the rotation angle θ.

rotation3D(θ, n=[0,0,1])

Returns SO(3) rotation matrix parametrized by rotation angle θ and rotation axis n=[n1,n2,n3].

Note: This routine will always normalize n.

signedangle(e1::T,e2::T; z=nothing)

returns the signed angle between 3D vectors e1 and e2. z should be the normal vector.

Lattice

Type that contains all information about a lattice.

Fields

• basis::Matrix: columns are basis vectors
• latticedim::Int: how many of the basis vectors are lattice vectors
• spacecoordinates::Matrix: columns are coordinates of orbitals in the unit cell (w.r.t. basis vectors)
• extracoordinates::Matrix: additional non-spatial coordinates (e.g., could be sublattice index)
• extralabels::Dict: labels for the non-spatial coordinates
• specialpoints::LabeledPoints: High-symmetry points of the lattice (see ?Structure.Paths.LabeledPoints)

Constructors

Lattice(basis::Matrix)
Lattice(basis::Matrix [, latticedim::Int], spacecoordinates::Matrix [, extracoordinates::Matrix]; extralabels=Vector{String}(), specialpoints=LabeledPoints())

Property functions

latticedim(lat::Lattice)
countorbitals(lat:Lattice)
spacedim(lat::Lattice)
extraspacedim(lat::Lattice)

hasdimension(lat::Lattice, name::String)
assertdimension(lat::Lattice, name::String)

basis(lat::Lattice, ...)
getA(lat::Lattice, ...)
getB(lat::Lattice, ...)

coordinates(lat::Lattice, ...)
positions(lat::Lattice, ...)
allpositions(lat::Lattice, ...)
extracoordinates(lat::Lattice, ...)

filterindices(lat::Lattice, name::String, condition::Function)

Method functions

addbasis!
addorbital!
addorbitals!
addextra!
setextracoordinates!
fractionalize!
foldfractional
foldcoordinates!
rotatebasis!
rotatecoordinates!
translate!
displace!
displaceZ!
mirrorZ!
mirrorZ
newdimension!
mergelattices!
mergelattices

_lattice_reduce(A) → (A_red, U)

Lagrange (Gauss) reduction of the lattice basis A (columns are lattice vectors): returns a reduced basis A_red = A * U spanning the same lattice, with U unimodular (integer, |det U| = 1). For a reduced 2D basis ‖b1‖ ≤ ‖b2‖ and |b1·b2| ≤ ½‖b1‖² (angle in [60°,120°]), so fixed-size integer search boxes are sufficient to enumerate near-origin lattice points — which is false for sheared inputs like [1 6; 0 1]. Only size(A,2) == 2 is reduced; other dimensions return U = I unchanged.

bistack(lat, δz; fracshift=[0.0, 0.0])

Create a bilayer by duplicating lat and offsetting the upper copy by δz along the third (z) coordinate. Optionally shift the second layer in fractional in‑plane coordinates by fracshift before stacking. Adds a "z" extra coordinate if missing. Returns the new stacked Lattice.

crop2unitcell!(lat)

Remove orbitals whose fractional coordinates fall outside [0,1) along the primitive directions. Operates in place and returns lat.

crop2unitcell(positions, Λ)

Return the subset of positions (Cartesian) that lie inside the primitive unit cell spanned by Λ (direct basis). Returns a matrix whose columns are the kept positions.

displace!(lat, f::Function)

Function f takes orbital i at position pi and displaces it by vector vi = f(p_i).

displaceZ!(lat, f::Function)

Function f takes orbital i at position pi and displaces it in the third coodinate by zi = f(p_i).

extralabelsinorder(lat) → Vector{String}

Extra-coordinate labels ordered by their row index in extracoordinates(lat). Use this (never collect(keys(lat.extralabels)), which is Dict hash order) whenever labels are paired with extracoordinate rows — e.g. when rebuilding a Lattice from an existing one.

fillregion(lat::Lattice, f::Function)

Takes the d-dimensional lattice and tiles the d-dimensional region defined by function f. Returns a 0-dimensional lattice (single unit cell with no periodicities).

f(p)=true if point p belongs to the region and f(p)=false otherwise

Note: f must define a finite region that includes the origin.

foldBZ!(points, lat::Lattice; shift=0.0)

Fold k‑points (columns of points) into the first Brillouin zone of lattice lat. K‑points are assumed to be in fractional (reciprocal‑basis) coordinates. An optional shift vector can be subtracted before folding.

foldPC!(points, lat::Lattice; shift=0.0)

Fold lattice coordinates (columns of points) into the first primitive unit cell of lat. Coordinates must be fractional in the direct‑lattice basis. An optional shift vector can be subtracted before folding.

foldPC!(lat::Lattice; shift=0.0)

In‑place variant that folds lat.spacecoordinates into the first primitive unit cell of lat. Coordinates are interpreted as fractional. Returns the folded coordinate submatrix view for convenience.

foldcell!(points::AbstractMatrix{Float64}, basis::Matrix{Float64})

Fold columns of points into the first Wigner–Seitz cell defined by the column vectors of basis (direct-lattice basis), then calls foldcell_fromneighbors! with the metric basis' * basis.

The half-plane set is built from all {-1,0,1} combinations of the Lagrange-reduced basis (mapped back to integer coefficients of basis), which contains every Voronoi-relevant vector of a 2D lattice. Using only the single shortest neighbor shell (the old behavior) never folds along directions longer than the shortest lattice vector — e.g. on a rectangular diag(1,2) cell, points were left unfolded along the long axis.

foldcell_fromneighbors!(points, gvectors[, M])

Fold columns of points into the first Wigner–Seitz cell using a list of neighbor-cell vectors gvectors (both in fractional coordinates). The optional metric M = B'B can be supplied to define distances in a non-orthogonal basis, where B = [G1 G2 …] collects the reciprocal lattice vectors as columns.

Arguments

• points::AbstractMatrix{Float64}: columns are coordinates to be folded (fractional units, i.e. in the basis of the lattice vectors).
• gvectors::Matrix{Float64}: columns are the neighbor-cell shift vectors that span the Wigner–Seitz region (fractional units).
• M::AbstractMatrix (optional): metric used to compute distances; defaults to the identity (orthonormal basis).

Returns the same points matrix with all columns folded in place. This is the low‑level implementation used by higher‑level foldcell! methods.

Calculate the dual lattice of lat.A. Here we use the general formula $B = A * (A^T * A)^-1$. That works also when the d-dim lattice is embedded in D-dim space.

getneighborBZ(lat, k=1; halfspace=true, innerpoints=false, excludeorigin=true)

This is the analogue of method getneighborcells(), except that it looks for nearest neighbor cells in reciprocal space.

getneighborcells(A, k=1; halfspace=true, innerpoints=false, excludeorigin=true)

Find the k-th-nearest neighboring unit cells given lattice vectors A[:,i], returned as integer coefficient vectors in the basis A. If halfspace=true, the list only contain [I,J] without its partner [-I,-J]. If innerpoints=true, returns all neighboring cells up to and including the k-th ones.

The enumeration runs over a fixed box in the Lagrange-reduced basis (see _lattice_reduce) and maps the results back, so skewed/non-reduced bases (e.g. sheared supercells) find their true nearest cells; a fixed box in the raw basis silently missed them.

getneighborcells(lat, k=1; halfspace=true, innerpoints=false, excludeorigin=true)

A naive implementation to find a list of k-th-nearest neighboring unit cells. If halfspace=true, the list only contain [I,J] without its partner [-I,-J]. If innerpoints=true, returns all neighboring cells up to and including the k-th ones.

moire_supercell(n, m=1)

Return the integer 2×2 supercell matrix S whose columns are the moiré Bravais vectors expressed in the original lattice basis: [n -n-m; n+m 2n+m], with abs(det(S)) = 3n²+3nm+m² primitive cells per layer.

This is a √3-rotated supercell relative to the smaller "minimal" moiré cell (area n²+nm+m²) commonly seen in TBG papers; LatticeQM uses the larger cell because the construction algorithm (in-plane parity flip + tile + rotate + crop, in twist) is only commensurate with the rotated layer when there is this extra √3 of room. A direct minimal-cell construction would need a different algorithm (no parity flip, direct integer-lattice enumeration) — not currently implemented.

repeat(spacecoordinates[, Λ], repeat=[0:0, 0:0])

Tile a set of point coordinates by integer translations of the lattice. Returns the concatenated coordinate matrix.

Arguments

• spacecoordinates::AbstractMatrix: columns are point positions in Cartesian space.
• Λ::AbstractMatrix (optional): direct lattice basis (defaults to 2D identity).
• repeat::Vector{<:AbstractRange}: ranges along a1 and a2, e.g. [0:1, 0:1] for 2×2.

Example

coords2x2 = repeat(coords, A, [0:1, 0:1])

sortposition!(lat, name::String[, sortfunc])

Sort orbitals by an extra coordinate name (e.g. "z", "sublattice"). Useful for layered plots so that layers appear in front/back order. Returns the permutation applied.

sortposition!(lat, index::Int[, sortfunc])

Index‑based variant of sortposition!, using the index‑th extra coordinate. Returns the permutation applied.

supercellpoints(M::AbstractMatrix{Int}; offset::Float64=sqrt(eps()))

Given an integer lattice Z^D and an integer matrix M describing a (possibly non‑orthogonal) superlattice, return the integer points inside one supercell of the new lattice. The number of columns equals abs(det(M)).

The small offset avoids boundary ambiguities for points lying exactly on cell faces.

superlattice(lat::Lattice, S::Matrix{Int}, supercellints::Matrix{Int}; kwargs...)

Low‑level constructor where the integer coordinates inside one supercell are pre‑computed and provided as supercellints (columns). Not intended for general users.

superlattice(lat::Lattice, superperiods; kwargs...)

Construct a superlattice from lat.

• superperiods: either a vector of integers (scales each direct basis vector) or an integer matrix S whose columns define the linear transformation of the primitive basis, i.e. the new basis is A * S.

Copies atoms into the new supercell and returns a new Lattice with updated coordinates and specialpoints preserved. Extra coordinates (e.g. sublattice, layer, z) are replicated.

twist(lat, n; kwargs...)
twist(lat1, lat2, n; z=3.0, m=1, verbose=false)

Build a commensurate twisted bilayer from the 2D triangular lattice(s) lat1, lat2 (must share the same Bravais vectors). The twist angle is twistangle(n; m=m); the inter-layer separation along z is z.

Inputs are not mutated.

Keyword arguments

• z::Float64=3.0 — inter-layer distance.
• m::Int=1 — second moiré index. Must satisfy gcd(n, m) = 1.
• verbose::Bool=false — print the chosen angle and supercell.

twistangle(n; m=1, degrees=false)

Commensurate twist angle for the moiré index pair (n, m). Uses the convention cos α = (3n² + 3nm + m²/2) / (3n² + 3nm + m²). Set degrees=true to return the angle in degrees instead of radians.

graphene()

Convenience alias for honeycomb(1.42) using a typical C–C bond length in Ångström units.

honeycomb(a=1.0)

Two‑site honeycomb (graphene) lattice with lattice constant a (distance between nearest neighbors equals a). Adds a "sublattice" extra coordinate with values 0 (A) and 1 (B).

honeycomb_AA(a=1.0, z=3.0)

AA‑stacked honeycomb bilayer with interlayer distance z.

honeycomb_AB(a=1.0, z=3.0)

AB (Bernal) stacking for a honeycomb bilayer with interlayer distance z.

honeycomb_ABC(a=1.0, z=3.0)

ABC (rhombohedral) stacked tri‑layer honeycomb. Adds "sublattice" and "layer" extra coordinates.

honeycomb_BA(a=1.0, z=3.0)

BA stacking variant for a honeycomb bilayer with interlayer distance z.

honeycomb_bilayer(a=1.0, z=3.0; δ=[0.0, 0.0])

Bernal‑stacked honeycomb bilayer (AB) with interlayer distance z. Optional in‑plane shift δ (fractional coordinates) moves the top layer prior to stacking. Extra coordinate "sublattice" is provided.

honeycomb_twisted(N, a=1.0, z=3.0; fold=true, verbose=false)

Commensurate twisted bilayer graphene (TBG) with twist index N. If fold=true the structure is wrapped to the primitive cell. Uses the legacy √3-rotated moiré supercell (3n²+3n+1 cells/layer); see Lattices.moire_supercell for why the smaller minimal cell is currently unsupported.

honeycomb_twisted_ABAB(N, a=1.0, z=3.0; fold=true, verbose=false)

Four‑layer twisted stack with ABAB stacking within layers prior to twisting.

honeycomb_twisted_ABBA(N, a=1.0, z=3.0; fold=true, verbose=false)

Four‑layer twisted stack with ABBA stacking within layers prior to twisting.

smoothdisplaceZ!(lat, δz_even=0.055, δz_odd=0.0; sharp=1)

Apply a smooth out‑of‑plane displacement pattern to a hexagonal bilayer lat parameterized by even/odd amplitudes. The optional sharp parameter controls the smoothness (higher values approach a sign‑like modulation).

square(a=1.0)

Square Bravais lattice with lattice constant a. Returns a Lattice with one orbital per cell and standard special k‑points (Γ, X, Y, M).

triangular(a=1.0)

Triangular Bravais lattice with lattice constant a. Returns a Lattice with one orbital per cell and triangular‑lattice special points (Γ, K, M, …).

triangular_supercell(a=1.0)

A 3‑site triangular supercell useful for toy models (e.g. Kekulé). Provides a "sublattice" extra label with values 1,2,3.

triangular_twisted(N, a=1.0, z=3.0; fold=true)

Commensurate twisted bilayer built from triangular lattices using twist index N (moiré periodicity). If fold=true, folds positions back to the primitive cell. Sets a compact set of special k‑points.

LatticeQM

Library for tight-binding models defined on (periodic) lattices, providing convenient functions to build the operators and to obtain bands, expectation values, topological indices, linear response coefficients and mean-field solutions.

Submodules

• Structure
• TightBinding
• Spectrum
• Operators
• LinearResponse
• Meanfield

Any of these modules can be further explored, e.g., with ?TightBinding

Usage examples

See folder examples of the package.

Structure

Provides the struct Lattices.Lattice to define and manipulate lattices, and Paths.DiscretePath to deal with discretized paths.

Check out the submodules:

• Lattices
• Paths
• Geometries

Lattices

Provides the type Lattice and methods that act on this struct, see ?Structure.Lattices.Lattice.

Geometries

Provides predefined lattice objects (such as two-dimensional honeycomb lattice).

Example

import Structure.Geometries

lat = Geometries.honeycomb_twisted(11)
plot(lat, 3; supercell=[0:1,0:1])

Parallel

Single primitive for k-space iteration across executors.

The package historically reimplemented if multimode == :distributed && nprocs() > 1 ... inside every routine that loops over k-points. This module replaces that with one parameterised loop:

Parallel.kspace_foreach!(body!, ks, exec; scratch_factory, progress)

body!(scratch, j, k) runs once per k-index. scratch is a per-task buffer (thread-local or worker-local), produced once per chunk by scratch_factory(), not per k. This is the difference that lets dense bandmatrix on a 364-orbital TBG cell drop from ~970 MB allocations to ~85 MB.

Supported executors:

• SerialExec() — no parallelism.
• ThreadedExec(n=Threads.nthreads(); schedule=:dynamic) — spawns n tasks on the default thread pool, each consuming work via a shared Channel. Dynamic by default, so irregular per-k cost (e.g. sparse Krylov convergence varying near band touchings) is load-balanced.
• DistributedExec(n=nworkers()) — uses pmap over chunks; one scratch per worker.

configure_blas!(exec) pins BLAS to one thread per worker for non-serial executors. Skipping this is a 100× footgun on a multi-core box (each worker's default BLAS pool oversubscribes the cores). Always call it once before running the executor.

Naive implementation of combining the linear spaces of two hopping models.

DistanceWindowHopping(a, halfwidth, t0)

Callable hopping rule: returns t0 when the bond length is in (a - halfwidth, a + halfwidth), otherwise zero(t0).

The struct is callable with three signatures (matching gethops's expectations): h(r::Real) # scalar distance h(r1::AbstractVector, r2::Real=0.0) # uses norm(r1) (truncated to 3D) h(r1::AbstractVector, r2::AbstractVector) # uses norm(r1 - r2) (truncated to 3D)

t0 may be any number type; the return type is T = typeof(t0) regardless of the branch, so downstream gethops/hops! see a type-stable function. Every call to nearestneighbor! produces a value of the same concrete struct type (parameterised by T), so the gethops → hops! → hoppingmatrix! pipeline compiles once per T instead of once per closure.

addhaldane!(hops, lat, t2; ϕ=π/2, spinhalf=false, cellrange=1, mode=:none, zmode=:none)

This method is a somewhat inefficient way to compute the haldane hopping matrix. The only upside to it is that it uses methods that I already implemented and that it is fairly general.

addsublatticeimbalance!(hops, lat, Δ; kwargs...)

Add a sublattice‑staggered chemical potential (imbalance) of magnitude Δ to hops on lattice lat. The on‑site energy is shifted by Δ*(s - 1/2) with sublattice coordinate s (A: s=0, B: s=1), so positive Δ lowers sublattice A by Δ/2 and raises sublattice B by Δ/2.

No‑op for Δ≈0.

addvalley!(hops, lat, fz=x->sign(x[3]+1e-3); kwargs...)

Add a valley mass term to hops on lat. Customise the layer sign via fz.

addzeeman!(hops, lat, Mv::Function)

Add a site‑dependent Zeeman coupling defined by the vector field Mv(r) to hops on lattice lat. The operator acts in spin‑1/2 space via σ·M(r) at each site. Modifies hops in place and returns it.

addzeeman!(hops, lat, M::AbstractVector; format=:dense)

Add a uniform Zeeman field M = [Mx, My, Mz] (constant across the lattice). If the norm of M is (numerically) zero, this is a no‑op. Returns hops.

blockmatrix(mat, I, J, V)
blockmatrix!(mat, I, J, V)

Writes a sparse matrix from block matrices.

N (Int): Size of the resulting N-times-N square matrix I,J (Vector{Int}): Row/Column coordinates vecofmats (Vector{Matrix}): Vector of block matrices (must be equal sized in the current implementation)

diophantine_cherns(p, q, r; Nw, smargin=0) -> Vector{Tuple{Int,Int}}

Candidate (s, C) labels of a Hofstadter gap from the TKNN / Středa–Wannier Diophantine equation

r = q·s − p·C,

for flux Φ=p/q (coprime) with r magnetic subbands filled below the gap and Nw orbitals per primitive cell. C ∈ ℤ is the gap Chern number as reported by Spectrum.gapcherns / hofstadter_cherns (the sign convention those routines produce), and s ∈ ℤ the band-filling label; the integrated density per primitive cell is n = r/q = s − C·Φ (in the conventional Wannier-diagram form n = s + t·Φ the Středa slope is t = −C).

The equation fixes C ≡ −p⁻¹ r (mod q) only, so every branch inside the heuristic window −smargin ≤ s ≤ Nw + smargin is returned, ordered by |C|.

reducelatdim(hops, index::Int)

Drop a lattice dimension. For example, a two-dimensional lattice with lattice vectors a1, a2 can be turned into a 1D ribbon with lattice vector a1 by dropping all hoppings along a2 (index=2) or inteo a 1D with lattice vector a2 by dropping all hoppings along a1 (index=1).

Tip: Together with superlattice(hops, periods) one can control the width of the finite ribbon before applying reducelatdim.

getcurrentoperators(lat, t::Function; kwargs...)

Build current operators directly from a hopping function t(r1, r2; ...). Returns one Hops per lattice direction. kwargs are forwarded both to t and to the Hops(lat, ...) constructor (e.g. cellrange, format).

getcurrentoperators(lat, hops::Hops)

Return the current operators [J_1, ..., J_d] (one per lattice direction) for the tight-binding Hamiltonian represented by hops on lattice lat.

Each J_α is an Hops whose (i,j) element at offset R is -i (r_i + A·R - r_j)_α · h_{ij}(R), i.e. ∂/∂k_α of the Bloch Hamiltonian with the standard position-included gauge.

gethops(lat::Lattice, t::Function; cellrange=1, format=:auto, vectorized=false)

Iterates over pairs of orbitals/atom positions (r1,r2) in lattice lat and evaluates the hopping elements t(r1+R,r2) for each lattice vector R.

By default, vectorized=false. For huge systems use vectorized=true and make sure the hopping function t accepts matrices as arguments. The keyword argument format can be :dense or :sparse. For :auto, small systems will be dense and huge problems are assumed to be sparse.

Returns the hopping elements in the format Dict(R => t_R)

getoperator(lat::Lattice, name::String, args...; kwargs...)

Return common projection/observable operators defined on lat by name.

Supported names (case-insensitive where noted):

• "sx", "sy", "sz" (or "MX/MY/MZ", "SX/SY/SZ"): Spin Pauli components.
• "spinup", "spindown": Projectors onto ↑/↓ in the local spin basis.
• "Sn": Spin along custom axis via Sn(lat, n̂) (see keyword forms).
• "layer": Layer projector for multilayer geometries.
• "sublattice": A/B sublattice projector; also "sublatticeA", "sublatticeB".
• "sublatticeAspin", "sublatticeBspin": Sublattice-resolved spin projectors.

For interaction kernels (e.g. "Hubbard", "CappedYukawa"), use the dedicated constructors in Operators/Meanfield (see gethubbard, getcappedyukawa).

Returns a Hops-like operator suitable for Spectrum.getbands projectors or expectation-value routines.

getzeeman(lat, M; kwargs...)

Construct a Zeeman term as a Hops operator. M can be either a constant magnetic field vector [Mx, My, Mz] (units absorbed in gμB) or a function Mv(r) returning a 3‑vector at position r for spatially varying fields.

Use addzeeman!(...) to add in place to an existing operator.

hofstadter(hops, lat, Q)

Determines the energie spectrum as function of rational magnetic flux $\Phi=p/q$, where p, q are coprime integers with 1<= q <= Q and 1<=p<q. Returns a list of fluxes and a list of energies at each flux.

In this implementation we evaluate at the $\Gamma$-point of the magnetic cell.

hofstadter_cherns(hops, lat, Q::Int, NX::Int=18; gaptol=1e-6, kwargs...)

Topological Hofstadter (Wannier) diagram: the Chern number of every spectral gap as a function of rational magnetic flux $\Phi=p/q$ per unit cell (coprime p, q with 1 ≤ q ≤ Q and 1 ≤ p < q).

For each flux the magnetic supercell Hamiltonian is built with peierlsoutplane and its gaps are labelled by Spectrum.gapcherns — the cumulative non-abelian Chern number of the magnetic subbands below each gap (= Hall conductance σxy in e²/h). The magnetic BZ is sampled on an NX × NX grid; NX ≳ 18 is recommended, as the densely packed subbands need a fine grid to converge.

Filling is reported per original unit cell, ν = n/q for n subbands filled, so the gaps fall on Diophantine/Středa lines ν = s − C·Φ with integer Chern number C and band offset s — the sign is fixed by the flux sign of peierlsoutplane together with the Berry orientation of Spectrum.gapcherns and is pinned by a regression test (honeycomb, all gaps at q ≤ 5).

Returns a NamedTuple of equal-length vectors (; flux, filling, energy, width, chern), one entry per (flux, gap): flux::Rational is Φ=p//q, filling is ν per unit cell, energy is the gap-centre energy, width is the gap size, and chern is the gap Chern number. With the default near-zero gaptol the list includes numerically tiny pseudo-gaps; filter or weight by width for a clean diagram.

kwargs are forwarded to the diagonalisation (multimode, executor, format).

hofstadter_dos(hops, lat, q_max::Int, (f_min, f_max), N=300)

See hofstadterdos(hops, lat, qmax, frequencies).

Returns a list of fluxes, frequencies and a the dos at each flux.

hofstadter_dos(hops, lat, q_max::Int, frequencies::AbstractVector; klin=100, Γ=0.05)

Determines the energie spectrum as function of rational magnetic flux $\Phi=p/q$, where p, q are coprime integers with 1<= q <= qmax and 1<=p<q. Note that qmax determines the size of the largest magnetic supercell.

The density of states (DOS) at each flux is calculated on a discrete k-grid (resolution give by klin). The parameter Γ determines the energy broadening when calculating DOS.

Returns a list of fluxes and a the dos at given frequencies at each flux.

hofstadter_gaplabels(hops, lat, Q::Int, NX::Int=18; gaptol=0.05, Nw=hopdim(hops), smargin=0)

Cheap (energies-only) Hofstadter gap labelling: the Diophantine/Středa Chern label of every spectral gap as a function of rational flux Φ=p/q per unit cell (coprime p, q with 1 ≤ q ≤ Q and 1 ≤ p < q).

For each flux the magnetic supercell is built with peierlsoutplane and its gaps located by Spectrum.spectralgaps (no eigenvectors, no Berry — much cheaper than hofstadter_cherns). Each gap is then labelled with diophantine_cherns: chern is the smallest-|C| candidate and nsol the number of branches inside the −smargin ≤ s ≤ Nw + smargin window. Nw is the number of orbitals per primitive cell. Every gap is reported: nsol == 0 means no in-window candidate exists (then chern is set to 0 and is meaningless — widen smargin), nsol > 1 flags a gap whose label needs disambiguation (a gap-family fit or a hofstadter_cherns anchor).

The labels are Level-1 guesses; see the diophantine_cherns warning for the systematic failure modes (out-of-window intercepts, Dirac Landau fans — both occur already for plain NN honeycomb). Treat chern as a candidate to be pinned by anchors, not as a result.

Filling is reported per primitive cell, ν = r/q. Returns a NamedTuple of equal-length vectors (; flux, filling, energy, width, q, r, chern, nsol).

This is the cheap Level-1 layer of the standard Hofstadter workflow; use hofstadter_cherns for the exact Berry-curvature anchors.

nearestneighbor!(hops, lat, t0=-1.0; a=1.0, halfwidth=0.01, kwargs...)

Add nearest-neighbour hopping t0 (between sites at distance ≈ a) to hops.

Backed by DistanceWindowHopping(a, halfwidth, t0); pass that struct directly to addhops! if you need the same hopping rule across multiple lattices without redundant compilation.

peierls!(hops, lat, B)

Add Peierls phases to operator hops on a finite (0-dimensional) lattice lat for the uniform magnetic field B=(B1,B2,B3) (in flux quanta per unit area), using the symmetric gauge; see uniformfieldphase for the sign convention.

A uniform field in a fixed gauge is not lattice-translation invariant, so this entry point refuses periodic lattices: the resulting Bloch Hamiltonian would be non-Hermitian and physically meaningless. For periodic systems use

• peierlsoutplane(hops, lat, p, q) — out-of-plane flux Φ=p/q with the correct magnetic supercell;
• peierlsinplane!(hops, lat, B) — in-plane fields (that gauge is translation covariant).

peierls!(hops, lat, phase)

This method implements Peierl's substitution on a tight-binding model by adding the appropriate phases to each hopping amplitude in the Hamiltonian given by hops.

The phase function phase must be a function with the signature phase(r1::AbstractVector, r2::AbstractVector)

Note:

• this should be the last step when constructing a tight-binding Hamiltonian
• You need to make sure that the phases that you add do not break translational symmetries. This is not a trivial matter and may lead to unexpected/undetected mistakes. If the phases of a bond and its conjugate partner do not pair up, the result is a non-Hermitian operator; this function warns (once) when that happens, because the dense diagonalization path silently projects onto the Hermitian part.

peierls!(hops, lat, B)

Add Peierls phases to operator hops on lattice geometry lat for the uniform in-plane magnetic field B=(B1,B2). Uses in-plane gauge, see uniformfieldphase_inplane(...).

peierlsoutplane(hops, lat, p, q)

Add uniform out-of-plane magnetic field B=(0,0,B3) to operator hops on lattice geometry lat, such that the flux per unit cell is Φ=p/q.

This method automatically constructs and returns the correct superoperartor and magnetic supercell.

setfilling!(H, kgrid::Structure.Mesh, filling; kwargs...)

Variant that accepts an explicit k-grid kgrid (e.g., from regulargrid). Computes μ on that grid and shifts H in place.

setfilling!(H, filling; nk=100, kwargs...)

Shift the chemical potential of Hamiltonian-like object H so that the resulting electronic filling matches filling (0–1 per spin). Internally uses Spectrum.chemicalpotential on a regular nk × nk grid and then applies addchemicalpotential!(H, -μ).

Keywords are forwarded to the underlying solver (e.g., temperature T, diagonalization format).

superlattice(hops, periods)

Turn a given hopping model into a superlattice model by copying cells and hoppings. This method can be useful as preparation before adding modulations that change the periodicity of the model.

uniformfieldphase(r1,r2; B)

The Peierls phase (in units of 2π) acquired by a hop from r2 to r1 in a uniform magnetic field B=(B1,B2,B3), computed in the symmetric gauge A = ½ B × r.

Convention

All Peierls routines in this module share the convention

t_ij → t_ij · exp(+i 2π/φ0 ∫_{r_j}^{r_i} A·dl),

i.e. this function returns +(1/φ0) ∫_{r2}^{r1} A·dl = ½ B·(r2 × r1). This matches uniformfieldphase_outplane (pinned by the Hofstadter/Chern tests), so the same nominal B produces the same field direction in every entry point. B is measured in flux quanta φ0=h/e per unit area and must be a 3-vector; 2D positions are treated as lying in the z=0 plane.

uniformfieldphase_inplane(r1,r2; B)

The Peierls phase (in units of 2π) acquired by a hop from r2 to r1 in a uniform in-plane magnetic field B=(B1,B2,0), in the gauge A = (B2 z, −B1 z, 0).

Returns +(1/φ0) ∫_{r2}^{r1} A·dl — the same sign convention as uniformfieldphase and uniformfieldphase_outplane; see the convention note in uniformfieldphase. This gauge depends only on the hop vector and the mean height z̄, so it is lattice-translation covariant and safe for periodic 2D lattices.

The magnetic field should be in units of flux quanta $\phi_0=h/e$ per unit area.

uniformfieldphase_outplane(lat, Φ)

Passes the lattice vectors a1 and a2 to uniformfieldphase_outplane(a1,a2,Φ).

uniformfieldphase_outplane(a1,a2)

This returns the phase function that respects translational symmetry along $a_1$ but not along $a_2$. For rational flux Φ=p/q, this gauge is periodic in a2'=q*a2.

The magnetic field should be in units of flux quanta $\phi_0=h/e$.

valley(lat; spinhalf=false, kwargs...)

Construct a valley operator on lat. If spinhalf=true, attach spin structure to match spinful Hamiltonians.

BandData

Struct to store BandData obtained from method getbands(...). The fields BandData.bands is a matrix of eigenvalues at each point along a path, the field BandData.obs is a matrix of expectation values and BandData.path is a DiscretePath object (contains discrete points and point labels).

Can be saved conveniently with savedlm(bands; path="data") and plotted with plot(bands).``

bandmatrix(H, ks::Matrix{Float} [, As]; kwargs...)

Calculates the energies for operator H(k) for each column vector k of matrix ks. If operators As=[A1, A2, ...] are given, their expectation values are calculated and stored for each eigenvector.

Accepts the same keywords as geteigvals, geteigvecs, geteigen. In particular: format (:sparse or :dense) and num_bands::Int.

Parallelism is selected via multimode:

• :serial — no parallelism
• :multithreaded / :threaded — Threads.@spawn over the default pool
• :distributed — pmap over Distributed.workers()
• :auto (default) — distributed if workers exist, else threads if >1, else serial

For full control, pass a Parallel.Executor directly via the executor kwarg (e.g. executor=Parallel.ThreadedExec(8; schedule=:static)).

Returns (bands, obs): bands is num_bands × num_k, obs is the matching tensor of projector expectation values.

Example

using LatticeQM
lat = Geometries.honeycomb()
h = Operators.graphene(lat)
ks = kpath(lat; num_points=200)
bands, obs = bandmatrix(h, ks.points, Operators.valley(lat))

berry(H, NX::Int, NY::Int=0, bandindices::AbstractArray=[])

Convenience method for berry(statesgrid).

berry(statesgrid)

Goes through each plaquette (i,j),(i+1,j),(i+1,j+1),(i,j+1) and calculates the (non-abelian) plaquette phase. statesgrid is a four-dimenional array, containing the discretization information and the occupied states.

You can create a statesgrid with statesgrid(H, nx, ny, bandindices) or use the wrapper berry(H, nx, ny, bandindices).

berryalongpath(H, kpoints)

Calculate the abelian Berry curvature for each band along a path of discrete k points. It builds little plaquettes along the path between the kpoints[:,i] and kpoints[:,i+1].

fermisurfacedensity(H, energies; broadening=:auto, lat=nothing, num_points=15, kwargs...)

High‑level convenience that computes band energies on a regular k‑grid and then returns (kgrid, density) as produced by fermisurfacedensity_fromdata. If lat is provided, the grid is folded into the first Brillouin zone and mapped to Cartesian reciprocal space via getB(lat).

fermisurfacedensity_fromdata(bands, energies, obs; broadening=:auto)

Weighted variant where obs provides an observable per band and k‑point of the same shape as bands. Each Lorentzian contribution is multiplied by the corresponding weight.

fermisurfacedensity_fromdata(bands, energies; broadening=:auto)

Broadened Fermi‑surface density without weights/observables. See fermisurfacedensity_fromdata(bands, energy; ...) for details.

fermisurfacedensity_fromdata(bands, energy; broadening=:auto)

Compute a broadened Fermi‑surface “density” from a band‑energy matrix bands of size (nbands, Nk) evaluated on a k‑grid. energy can be a scalar or a vector of Fermi energies. The result has size (length(energies), Nk) and sums contributions from all bands using a Lorentzian of width broadening.

If broadening == :auto (default), a heuristic width based on the band dispersion is used.

gapcherns(H, NX::Int, NY::Int=NX; gaptol=1e-6, multimode=:auto, executor=nothing, format=:dense)

Chern numbers of the spectral gaps of a two-dimensional Bloch Hamiltonian H.

The gaps are located with spectralgaps; for each gap the returned Chern number is the cumulative, non-abelian Berry flux of the occupied manifold (the lowest n bands),

C(n) = sum(berry(states[:, :, :, 1:n])),

i.e. the Hall conductance σxy (in units of e²/h) when the chemical potential lies in the gap. This manifold quantity is gauge-robust even where individual bands touch — unlike the single-band Chern numbers from getcherns, which are ill-defined for the touching subbands of e.g. a Hofstadter spectrum.

Returns a Vector of NamedTuples (; n, elo, ehi, chern), ordered by energy: n bands lie below a gap spanning [elo, ehi] with Chern number chern.

The grid must resolve the Berry curvature; for the densely packed magnetic subbands of a Hofstadter problem NX ≳ 18 is a safe default (coarser grids can misreport wide-Chern narrow gaps).

See also spectralgaps, getcherns, berry, Operators.hofstadter_cherns.

getbands(H, ks::Union{DiscretePath, AbstractMatrix} [, As]; kwargs...)

Calculates the bands for operator H along discrete path ks and if operators As=[A1, A2, ...] are given, their expectaction values are calculated and stored for each eigenvector.

Note that ks is a discrete path object as returned by kpath(lat::Lattice,...).

Accepts the same keywords as geteigvals, geteigvecs, geteigen. In particular: format (:sparse or :dense) and num_bands::Int.

Returns a BandData object (with fields bands, obs, path).

Example

using LatticeQM

lat = Geometries.honeycomb()
h = Operators.graphene(lat)
ks = kpath(lat; num_points=200)
valley = Operators.valleyoperator(lat)

bands = getbands(h, ks, valley)

using Plots
plot(bands)

getberry!(bands, h, ks)

Calculates and appends the data from berryalongpath(h,ks) to the data object bands. This is a convenience method that is useful when plotting band diagrams with Berry curvatures colored it.

getbroadening(broadening, bands)

Internal helper that resolves broadening to a numeric value. For :auto, uses the mean band standard deviation divided by √Nk.

getcherns(wavefunctions::Function, NX::Int, NY::Int=0, bands::AbstractArray=[])

Returns the chern numbers of wavefunctions corresponding to all bands in bands, where NX and NY denote the coarseness of discretization of k-space. For bands=[] it returns all chern numbers.

getdos(h, emin::Float64, emax::Float64, num=500; kwargs...)

Computes the density of states of operator h(k) on the entire Brillouine zone, discretized on a grid with $ k{lin} \times k{lin} $ points. and for the frequencies ωs=(ωmin, ω2, ..., ωmax). The paremter $\Gamma$ is the energy broadening.

Accepts the same kwargs as getdos(h, ωs; klin, Γ, kwargs...).

Note: the current implementation only works for a two-dimensional Brillouine zone. Might change in the future, but for now use dos(h, ks, ω; Γ) syntax if needed.

getdos(h, ks, ω; Γ, parallel=true, format=:auto)

Computes the density of states of operator h(k) using the points ks=(k1,k2,...) and for the frequencies ω=(ω1, ω2, ...). The paremter $\Gamma$ is the energy broadening.

Mode can be :distributed or :serial, format can be :auto, :sparse or :dense.

getdos(h, ωs; klin, Γ, kwargs...)

Computes the density of states of operator h(k) on the entire Brillouine zone, discretized on a grid with $ k{lin} \times k{lin} $ points. and for the frequencies ωs=(ω1, ω2, ...). The paremter $\Gamma$ is the energy broadening.

Accepts the same kwargs as dos(h, ks, ω).

Note: the current implementation only works for a two-dimensional Brillouine zone. Might change in the future, but for now use dos(h, ks, ω; Γ) syntax if needed.

getwindnum(wavefunctions::Function, NX::Int, NY::Int=0, bandnr::Integer=1)

Returns the winding number (according to Rudner et al. 2013) of the band which is on position bandnr in the spectrum. NX and NY denote the coarsness of the discretization in k-space and wavefunctions is a function returning the wavefunctions of the problem. (NB: The winding number of a band is the sum of the chern numbers of all the bands below it including its own chern number. For undriven systems or systems which are driven but "normal" it is equal to the chern number of said band. It only makes sense when looking at Floquet bands ie quasi-energies.)

ldos(H, ks, frequencies; Γ=0.1, format=:sparse, kwargs...)

Site-resolved local density of states, k-averaged over ks and summed over the frequency window frequencies.

plaquettephase(S00, S10, S01, S11)

Calculates the phase the (non-abelian) phase winding around a plaquette.

This method is not meant to be called directly, it is used by berry(statesgrid).

spectralgaps(H, NX::Int, NY::Int=NX; gaptol=1e-6, multimode=:auto, executor=nothing, format=:dense)

Global spectral gaps of a two-dimensional Bloch Hamiltonian H, found on an NX × NY k-grid from eigenvalues only (no eigenvectors) — the cheap companion to gapcherns.

A gap below band n is recorded when the grid-minimum of band n+1 exceeds the grid-maximum of band n by more than gaptol (in the energy units of H).

Returns a Vector of NamedTuples (; n, elo, ehi), ordered by energy: n bands lie below a gap spanning [elo, ehi].

statesgrid(H, NX::Int, NY::Int=0, bandindices::AbstractArray=[]; kwargs...)

Evaluates the eigenvectors on a discretized grid (2D Hamiltonian only!) and stores the result (preserving the grid information). Useful when plaquette phases need to be calculated.

statesgrid[i,j,:,k] is the k-th eigenvector at gridpoint i,j.

Accepts multimode/executor and format like bandmatrix; defaults to :auto and :dense.

opticalconductivity(frequencies, H, J1, J2, ks; μ=0.0, Γ=0.025, T=0.1, ...)

Evaluate the Kubo current-current correlator at the set of frequencies and k‑points ks (columns). H, J1, and J2 are callable with a k‑vector and return the Hamiltonian and current operators, respectively.

Returns a complex vector χ(ω) normalized by the number of k‑points N:

χ(ω) = -(i/N) Σ_k Σ_{m≠n} (f_m - f_n) ⟨m|J1|n⟩⟨n|J2†|m⟩ / (ϵ_m - ϵ_n - ω - iΓ)

where f = fermidirac(ϵ; μ, T) and |m⟩, ϵ_m are the Bloch eigenstates at k. This is not yet the optical conductivity. Obtain σ(ω) by subtracting the static (diamagnetic) piece and dividing by the photon energy:

σ(ω) = -(χ(ω) - χ(0)) / (ω + iΓ)

as done in test/test_linearresponse.jl and extra/examples/graphene/opticalconductivity.jl.

μ is the chemical potential, Γ the phenomenological broadening and T the temperature (all in the same energy units as H). Additional keyword arguments are forwarded to the eigenvalue solver.

opticalconductivity(frequencies, H, J1, J2; klin, kwargs...)

Compute the Kubo current-current correlator χ(ω) on a regular klin × klin k‑grid using currents J1, J2. Convenience front‑end to the explicit‑grid method, whose docstring gives the formula and the χ(ω) → σ(ω) conversion.

opticalconductivity(frequencies, i, j, H, lat; kwargs...)

Compute the Kubo current-current correlator χ_ij(ω) for Hamiltonian H(k) using the Cartesian components i, j of the current operator defined on lat. Internally obtains J = getcurrentoperators(lat, H) and forwards to the low‑level routine, whose docstring gives the formula and the conversion from the returned χ(ω) to the optical conductivity σ(ω).

Common keywords: μ (chemical potential), Γ (broadening), T (temperature), plus diagonalization options forwarded to the eigen solver.

opticalconductivityXX(frequencies, H, lat; kwargs...)

Convenience wrapper computing the xx component of the Kubo current-current correlator χ_xx(ω). Builds current operators from H and lat and evaluates at the provided frequencies.

Note: despite the name, this returns the raw correlator χ(ω), not the optical conductivity σ(ω). See the low-level opticalconductivity(frequencies, H, J1, J2, ks; ...) docstring for the formula and the χ → σ conversion.

opticalconductivityXY(frequencies, H, lat; kwargs...)

Convenience wrapper computing the xy component of the Kubo current-current correlator χxy(ω) — not σxy(ω) itself; see the low-level opticalconductivity(frequencies, H, J1, J2, ks; ...) docstring for the χ → σ conversion. See opticalconductivity(frequencies, i, j, H, lat; ...) for details.

HartreeFock(h, v, μ=0.0; hartree=true, fock=true)

Mean-field functional for density (Hartree) and exchange (Fock) channels built from a base Hamiltonian h and interaction kernels v. Calling the struct on ρ updates the effective mean-field operator hMF and the scalar energy fields ϵH, ϵF (see Energy decomposition below). The SCF driver fills in ϵband and ϵkin once a density-matrix iteration has run.

h and v may use different matrix backends (e.g. dense h with sparse Hubbard v); they only need to share the lattice key type and dimensions.

Convention

v and h (and the resulting density matrix ρ) must use the same orbital basis. In particular, if your model has spin, the basis must already include it (e.g. via TightBinding.addspin); the Hartree term reads diagonal occupations from ρ and assumes those are the densities v couples to. For a spinful Hubbard U n_↑ n_↓ model, build v so its matrix elements between spin-↑ and spin-↓ orbitals encode U, not the same-spin diagonal.

Energy decomposition (physical-sign convention)

After hf(ρ):

• ϵH = +½ nᵀ V₀ n = ½ Tr[V_H ρ] — Hartree energy (positive for repulsive V₀)
• ϵF = -½ Σ_L Σ_{ij} v_L[i,j] |ρ_L[i,j]|² = ½ Tr[V_F ρ] — Fock / exchange energy (negative for repulsive v)

After one SCF iteration the driver sets:

• ϵband = Tr[hMF · ρ] — sum of occupied band energies
• ϵkin = Tr[h · ρ] = ϵband - 2(ϵH + ϵF) — kinetic / bare-Hamiltonian expectation value (derived from the variational identity E_band = ⟨T⟩ + 2(E_H + E_F))

The total HF energy can be read off in three equivalent ways — useful as self-consistency checks at convergence:

• E_HF = ϵkin + ϵH + ϵF (physical, additive)
• E_HF = ϵband - ϵH - ϵF (band − double-counting)
• E_HF = ½ (ϵkin + ϵband) (½ Tr[ρ (h₀ + hMF)])

MeanfieldGenerator

Abstract supertype for mean‑field functionals that map a density matrix ρ to an effective single‑particle Hamiltonian and scalar energy contributions. Concrete implementations include HartreeFock and HartreeFockBDG.

_scf_driver!(ρ0, ρ1, mf, update_ρ!; kwargs...)

Common SCF skeleton shared by solveselfconsistent! (k-space diagonalization) and solveselfconsistent_purification! (real-space canonical purification). The caller supplies a closure update_ρ!(ρ1, mf) -> ϵ_band that overwrites ρ1 from the current mean-field operator and returns the band-energy contribution Tr[H_MF · ρ_{k+1}] at this iteration. The driver caches that on mf.ϵband and derives mf.ϵkin = ϵband - 2(ϵH + ϵF) (variational identity), so callers read the full physical decomposition (ϵkin, ϵH, ϵF, ϵband) from the returned mean-field generator.

Keywords:

• convergenceerror::Bool=false — throw if SCF did not converge.
• callback=(ρ1 -> nothing) — invoked with the freshly computed ρ1 at each iteration. Existing user hook.
• log_callback=nothing — if non-nothing, called every iteration. The driver invokes it as either log_callback(iter, ϵ, residual, info) (4 args) or the legacy log_callback(iter, ϵ, residual) (3 args), picked via applicable. info is a NamedTuple (ϵkin, ϵband, ϵH, ϵF, μ).
• residual_fn=nothing — closure (ρ1, ρ0) -> Real that overrides the convergence metric used by fixedpoint!. Defaults provided by the caller (typically the commutator norm ‖[H_MF, ρ]‖_F).
• iterations::Int=500, tol::Real=1e-7, verbose::Bool=false — passed through to fixedpoint!.
• Any other kwargs are forwarded to fixedpoint! (e.g. β, p_norm, relative, show_trace).

solveselfconsistent!(ρ0, ρ1, mf, filling, ks; kwargs...)
solveselfconsistent!(ρ0, mf, filling, ks; kwargs...)
solveselfconsistent!(ρ0, ρ1, mf, filling; klin, kwargs...)
solveselfconsistent!(ρ0, mf, filling; klin, kwargs...)

Search for a self-consistent mean-field solution at the given filling (between 0 and 1). The MeanfieldGenerator mf (e.g. HartreeFock) maps a density matrix ρ to the effective single-particle Hamiltonian. k-space is discretized with the supplied points ks (or built from klin).

Returns (ρ, ϵ_GS, mf, converged, residual):

1. converged density matrix,
2. mean-field ground-state energy,
3. the updated MeanfieldGenerator (with chemical potential mf.μ),
4. convergence flag,
5. final residual.

Keywords (all optional): iterations, tol, T (electronic temperature), β (linear mixing), convergenceerror, multimode, verbose, hidebar, callback, log_callback. See _scf_driver! and fixedpoint! for details.

parallel=true (via multimode) helps when per-k diagonalization dominates (e.g. twisted bilayer graphene). For small problems, the serial path can be faster due to communication overhead.

solveselfconsistent_purification!(ρ0, ρ1, mf, filling, ks; kwargs...)

Real-space canonical-purification analogue of solveselfconsistent. Each SCF step replaces k-space diagonalization with a fixed-step McWeeny-style purification of the current mean-field Hamiltonian. Useful for very large systems where per-k diagonalization is the bottleneck.

Returns (ρ, ϵ_GS, mf, converged, residual). Same kwargs as solveselfconsistent!.

commutator_kspace_norm(H, ρ, ks; multimode=:auto, kweights=nothing) → Real

Frobenius norm of the commutator [H, ρ], evaluated in k-space:

\[‖[H, ρ]‖_F = \sqrt{\sum_k w_k\,‖H(k)\,ρ(k) - ρ(k)\,H(k)‖_F^2}\]

H and ρ must share an orbital basis (same matrix dimension per Hops block); their key sets may differ. Used as the SCF convergence metric: self-consistency is exactly [H_MF[ρ], ρ] = 0, so the commutator norm is a direct, basis-independent measure of how far the iteration is from a fixed point.

doublecounting(hf) → Float64

Total interaction double-counting correction of the mean-field band energy: E_GS = ϵband − doublecounting(hf) and ⟨T⟩ = ϵband − 2·doublecounting(hf). For plain Hartree-Fock this is ϵH + ϵF; BdG generators add the pairing channel ϵP (the quasiparticle band energy contains ⟨½(c†Δc† + h.c.)⟩ = 2ϵP while the physical pairing interaction energy is ϵP).

enrichkeys!(ρ::Hops, klin::Int) → Hops

Expand keys(ρ) to cover the full Wigner-Seitz cell of the inverse k-grid that regulargrid(nk=klin^d) produces, so the SCF can populate every Fourier mode of ρ_full(k). Keys are added with zero matrix blocks; any existing block content is preserved.

This is the prerequisite for using residual_norm=:commutator as a hard convergence criterion: with sparse keys(ρ) (the default for e.g. on-site Hubbard v), bloch!(ρ, k) reconstructs only a low-rank approximation of ρ_full(k), and ‖[H, ρ]‖_F hits a non-vanishing truncation floor.

klin parity

klin should be odd. The WS cell {-Lmax, …, Lmax}^d for odd klin is closed under L → -L, so the lattice-Hops Hermiticity condition ρ[L]^† = ρ[-L] is satisfied without adding redundant keys. With even klin, hermitianize! adds keys at the BZ boundary that double-count the k = ±½ Fourier mode, leaving a small (~10⁻²) residual floor.

Cost

Memory: klin^d × N² complex doubles. For klin = 9, d = 2, N = 4 (graphene-Hubbard): 81 × 16 ≈ 1300 doubles. Negligible. For klin = 9, N = 1588 (TBG-N=11): ~ 200 MB — non-trivial but tractable.

Example

ρ_init = Meanfield.initialguess(v, :antiferro; lat=lat)
Meanfield.enrichkeys!(ρ_init, 9)  # commutator residual now valid
ρ, ϵ, mf, conv, res = Meanfield.solvehartreefock(
    hops, v, ρ_init, 0.5; klin=9, residual_norm=:commutator, tol=1e-7)
# converges in ~8 iterations to res ≈ 1e-9

fixedpoint!(f!, x1, x0; iterations=500, tol=1e-7, β=1.0, p_norm::Real=2,
            relative=true, acceleration=:linear, anderson_depth=5,
            show_trace=false, clear_trace=false, verbose=true,
            log_callback=nothing, residual_fn=nothing)

Damped fixed-point iteration with optional Anderson (Type-II) acceleration. f!(x1, x0) should overwrite x1 with the next iterate computed from x0 and may return a scalar (e.g. ground-state energy at this step).

After each call to f!:

• Linear mixing (acceleration=:linear, default): the new iterate is β * f(x_k) + (1 - β) * x_k. β=1 is undamped.
• Anderson acceleration (acceleration=:anderson, opt-in): keeps the last anderson_depth (x_i, f(x_i)) pairs and solves a small least-squares problem to combine them. β acts as a damping parameter (β=1 is undamped Anderson). Falls back to plain mixing on the first step (no history yet) and whenever the LS solve produces a non-finite update. Typical SCF problems converge ~5–10× faster than linear mixing on well-conditioned interacting systems. Caveat: at U ≈ 0 and T = 0 the residuals between iterations are nearly collinear, the LS becomes ill-conditioned, and γ blows up — Anderson is opt-in rather than the default for that reason. A future version may add automatic rank-revealing QR pruning to make it safe by default.

The default residual is ‖x1 - x0‖ / max(‖x1‖, eps()) if relative=true, else absolute. Pass residual_fn=(x1, x0) -> Real to override it — the callback is evaluated immediately after f!, before mixing rewrites x0, so it sees the raw (x_{k+1}^{f}, x_k) pair. The Anderson least-squares math always uses the iterate-difference form internally regardless of residual_fn. Set log_callback=(iter, ϵ, residual) -> ... to log every step without enabling the progress bar.

Sanity-checked against f_a(x) = 1/2 * (a/x + x), fixed point √a.

fock_energy(hf, ρ) → Real

Physical Fock (exchange) energy -½ Σ_L Σ_{ij} v_L[i,j] |ρ_L[i,j]|² (negative for repulsive v).

hartree_energy(hf, ρ) → Real

Physical Hartree energy +½ nᵀ V₀ n (positive for repulsive V₀).

mapspindensitymatrix(vs::AbstractVector, Is::AbstractVector{Int}, N::Int) --> SparseMatrix{Complex}
mapspindensitymatrix(vs::AbstractVector, Is::AbstractVector{Int}) = mapspindensitymatrix(vs, Is, length(Is)) --> SparseMatrix{Complex}
mapspindensitymatrix(vs::AbstractVector) --> SparseMatrix{Complex}

From a list of spin orientations vs (and optionally lattice indices Is, and optionally total number of lattice sites N), a spin density matrix for multiple lattice sites is generated.

solvehartreefock(h, v, ρ_init, filling; kwargs...)

Convenience wrapper around solveselfconsistent that constructs a HartreeFock functional from base Hamiltonian h and interaction kernel v. Returns the converged mean‑field solution and metadata.

solveselfconsistent(ρ0, mf::MeanfieldGenerator, filling, ks; kwargs...)
solveselfconsistent(ρ0, mf::MeanfieldGenerator, filling; klin, kwargs...)

Non‑mutating convenience wrappers around solveselfconsistent! that copy the initial density matrix ρ0, iterate the mean‑field functional mf (e.g., HartreeFock), and return the converged result together with energy and state.

The filling sets the target electronic filling (0–1 per spin). Supply either an explicit k‑grid ks or a grid resolution via klin (uses klin×klin).

Common keywords: iterations, tol, T (temperature), β (mixing), multimode (parallel), log_callback, residual_norm.

Returned energy decomposition

The returned mf carries the converged decomposition (physical-sign convention):

• mf.ϵkin — kinetic energy Tr[h₀ · ρ] (derived from the variational identity ϵkin = ϵband - 2(ϵH + ϵF), exact at the SCF fixed point)
• mf.ϵH — Hartree energy +½ nᵀ V₀ n (positive for repulsive V₀)
• mf.ϵF — Fock energy -½ Σ_L Σ_{ij} v_L[i,j] |ρ_L[i,j]|² (negative for repulsive v)
• mf.ϵband — band energy Tr[H_MF · ρ] = sum of occupied eigenvalues

The total HF energy (returned as the second tuple element) can be expressed in three equivalent ways at the SCF fixed point:

• E_HF = ϵkin + ϵH + ϵF (additive physical decomposition)
• E_HF = ϵband - ϵH - ϵF (band − double-counting)
• E_HF = ½ (ϵkin + ϵband) (½ Tr[ρ (h₀ + H_MF)])

When log_callback is supplied, the driver invokes it as log_callback(iter, ϵ, residual, info) if the callback can take a 4th arg (detected via applicable), with info = (ϵkin, ϵband, ϵH, ϵF, μ). Otherwise it falls back to the legacy 3-arg form log_callback(iter, ϵ, residual).

spindensitymatrix(d::Vector=[1,0,0]) --> Matrix{Float64}
spindensitymatrix(x::Real=1) --> Matrix{Float64}
spindensitymatrix(s::Symbol) --> Matrix{Float64}

Returns a spin-1/2 density matrix for given spin orientiation d=[dx,dy,dz]. Instead, one can also specify the configuration symbolically, in which case s should be one of :up,:down,:upx,:downx,:upy,:downy.

If vector has norm larger one it gets normalized, otherwise it is left unscaled.

spinspiraldensitymatrix(lat::Lattice, superperiods; n::Vector=[0,0,1], v0::Vector=[1,0,0])

Create a spin spiral density matrix with periodicity given by a superlattice characterized by superperiods' (see Structure.superlattice for details). The rotation goes around the fixed axisn=[n1,n2,n3]. The vectorv0` specifies the spin direction in one site.

Useful to create initial states for mean field calculations.

BdGOperator(h::Hops)
BdGOperator(h::Hops, Δ::Hops)

Bogoliubov–de Gennes operator wrapper that lifts a normal-state hopping object h to Nambu space and optionally adds a pairing block Δ. The resulting operator behaves like a Hamiltonian H(k) and can be passed to spectrum and mean-field routines.

The Nambu Hamiltonian built is

H_BdG(k) = [ h(k) − μ I    Δ(k)         ]
           [ Δ†(k)         −h(−k)^T + μ I ]

with the spinless / pre-spin-doubled convention Ψ_i = (c_i, c_i^†)^T. The hole block is constructed as −conj(h(R)), which equals −h(−R)^T only when h is Hermitian — the constructor asserts this. The chemical potential is added later via addchemicalpotential!(H::BdGOperator, μ), not at construction time.

Notes on the dimensions:

• hopdim(H::BdGOperator) returns the electron-sector size N (the number of orbitals in the underlying h).
• Spectrum.dim(H::BdGOperator) returns the Nambu size 2N (what the eigensolver sees).

Related helpers:

• addpairing!(H::BdGOperator, Δ) — insert/update the pairing block.
• addelectronsector!(H::BdGOperator, h) — add to the electron block.
• getelectronview / getpairingview — views into electron/pairing sectors.

getdensitymatrix!(ρ::BdGOperator, H::BdGOperator, ks, μ=0; kwargs...)

Compute the full Nambu density matrix ρ from the BdG Hamiltonian H at chemical potential μ. Internally shifts H by −μ while running the sum, then restores H to its prior state (try/finally-protected).

Returns the physical grand-canonical kinetic energy ⟨H_phys − μ N̂⟩, recovered from the BdG sum via

⟨H_phys − μ N̂⟩ = (1/2) Tr(ρ_BdG · H_BdG^shifted) + (Tr h_e − N μ) / 2,

where the constant term is the Nambu c-number from anticommuting c c† = 1 − c†c in the hole sector. Without it the reported value would equal twice the kinetic plus a μ- and N-dependent offset; this matters whenever ϵ_GS = ϵ_kin + ϵ_MF is consumed (e.g. comparing total energies between SCF runs at different parameters).

Spectrum.chemicalpotential(H::BdGOperator, ks, target::Real; T=0.01, ...)

Chemical potential at which the BdG Hamiltonian H reaches the target electron filling. Bisects on μ using Spectrum.filling(H::BdGOperator) above (i.e. the full BdG, including Δ). The default Spectrum.chemicalpotential(::Any, ...) strips Δ before fitting, which is fine for weak coupling but wrong for strong pairing — this method overrides that for BdGOperator so the SCF self-consistently solves μ even with a non-zero gap.

Spectrum.filling(H::BdGOperator, ks; μ=0.0, T=0.01, kwargs...)

Electron filling fraction at chemical potential μ for the BdG Hamiltonian H, computed from the full BdG density matrix (i.e. accounting for any pairing block Δ). This is Tr(ρ_el(k=0)) / N averaged over ks weights.

multimode defaults to :serial here for two reasons:

1. The BdG getdensitymatrix! mutates H in place via addchemicalpotential!(H, ±μ) to apply the chemical potential. The shift is bracketed around the k-loop (added before, undone after), so threaded reads are safe — but the bracket is exception-fragile if a worker throws. Serial keeps the contract simple.

2. Bisection in chemicalpotential(::BdGOperator) calls this in a tight loop (Roots.find_zero); the per-call threading overhead exceeds the per-call eigen cost on small (e.g. 144-k) grids. Pass multimode=:multithreaded or :distributed for larger workloads.

addpairing!(H::BdGOperator, Δ::Hops)

Add or update the pairing block of a Nambu Hamiltonian. The user must supply a Δ whose key set is closed under negation — every R in keys(Δ) must have a partner −R. The block built at offset R is

[ 0       Δ[ R] ]
[ Δ[−R]† 0     ]

Note that Δ[ R] and Δ[−R] parametrise the independent upper-right blocks at offsets R and −R (the lower-left at R is Δ[−R]†). The resulting BdG is Hermitian as a Hops regardless of any relationship between Δ[ R] and Δ[−R] — the only requirement is that both keys are present so the construction can run. This deliberately matches what the SCF produces: the mean-field pairing ΔMF[L] = v[L] .* ρΔ[L] is generically not dict-Hermitian (ρΔ is the upper-right block of ρ_BdG, which does not equal the lower-left in general), but the BdG built from it is still Hermitian.

periodicDrive(omega::Number=1)

Return an empty periodic drive with frequency omega.

mutable struct periodicDrive

Type representing a time-periodic potential applied to a lattice. Contains the fields:

• omega: Frequency of the potential.
• operators: An array of operators (Matrices) representing the Fourier components of the potential.
• ns: An array of integers containing the modes which are not zero. (Should not contain 0) This implies that operators and ns must have the same length and must be ordered such that operators[i] corresponds to the ns[i] fourier component for all i.

Example: The following returns a drive corresponding to a periodic potentisl V(t): V(t) = [0 1; -1 0]exp(-im2.5t) + [0 -1; 1 0]exp(3im2.5*t)

drive = periodicDrive(2.5, [[0 1; -1 0], [0 -1; 1 0]], [-1,3])

periodicDrive(omega::Number, funcFourier::AbstractArray{<:Number, 1}, operator::AbstractMatrix{<:Number})

Allows for initializing a periodicDrive using only one operator and Fourier components of a periodic function denoting the amplitude (over time) with which the operator is applied.

addFreq!(drive::periodicDrive, operator::AbstractMatrix, n::Integer)

Add a single harmonic mode to an existing periodicDrive.

addcos!(drive::periodicDrive, operator::AbstractMatrix, n::Integer)

Add the operator applied with a cosine of frequency n*omega to an existing periodicDrive.

getFloquetMatrix(M::Integer, H0::AbstractMatrix, drive::periodicDrive; hbar::Number=1.)

Build the effective Hamiltonian for several harmonic driving modes, i.e., $V(t) = \sum_j V_j \exp(i*j*n_j*\omega)$.

M: truncation H0: time averaged Hamiltonian (i.e., has already absorbed the constant m=0 mode), must be matrix drive: periodicDrive object

The drive V(t) is assumed Hermitian, which implies $V_{-n} = V_n^\dagger$ for its Fourier components. Harmonics whose partner -n is not stored in the drive are completed automatically with the adjoint, so the Floquet matrix is Hermitian even if only one member of a ±n pair is given. If both +n and -n are stored they must be adjoints of each other (checked within tolerance, otherwise an ArgumentError is thrown); the pair is then symmetrized so the result is exactly Hermitian.

getfirstFBZ(bands::AbstractArray, M::Integer, omega::Real) Takes a two dimensional array denoting a list of bands which came from diagonalizing a Floquet Hamiltonian and returns the bands lying in the first Floquet Brillouin zone, i.e., for each k point the dim quasienergies inside the window (-omega/2, omega/2], where dim is the dimension of the Bloch Hamiltonian. If the number of eigenvalues inside the window differs from dim (degeneracies at the zone edge or truncation artifacts), the dim states closest to the window are selected instead and folded back into (-omega/2, omega/2].

getfirstFBZ(bands::AbstractArray, M::Integer) Legacy method kept for backward compatibility: selects the middle dim rows of the energy-sorted bands. This is only correct as long as no quasienergy band crosses the edge of the Floquet Brillouin zone; otherwise a replica copy is included and the true band dropped. Pass the drive frequency omega (see the three-argument method) to get the correct windowed selection.

keepfirstFBZ!(data::Spectrum.BandData, H::FloquetOperator)

Throw away all band data except for the first Floquet Brillouin zone, i.e., the quasienergies inside (-omega/2, omega/2] with omega the drive frequency.

Converts an eigenvector vec obtained from solving the infinite (truncated at M) Floquet problem with drive frequency omega to an eigenvector in time space (at time t).

transformeigvecs(U::AbstractMatrix, M::Integer, omega::Real, t::Number) Converts a matrix of eigenvectors U obtained from solving the infinite (truncated at M) Floquet problem with drive frequency omega to eigenvectors in time space (at time t).

auto() -> Executor

Pick the most parallel executor that's actually available right now. Distributed wins if workers exist, else threads, else serial.

configure_blas!(exec; threads=nothing, verbose=true) -> Int

Pin BLAS threads for exec. Idempotent per target: the first explicit choice sticks (internal calls pass threads=nothing, which reuses the already-pinned value), workers are (re)pinned whenever their recorded target is missing or stale — including workers added after earlier calls.

• SerialExec → no change.

• DistributedExec → BLAS = threads (default 1) on master and every worker. Each Julia process is its own worker, so BLAS multi-threading oversubscribes.

• ThreadedExec → BLAS = threads (default 1). The library defaults to 1 because each Julia thread already calls LAPACK; with N Julia threads × M BLAS threads the box oversubscribes.

  This is workload- and backend-dependent. Empirically on Apple Silicon with the n=10 TBG cell (dense Hermitian heevr, 2648 orbitals): OpenBLAS prefers threads=3 (~30% better than threads=1); Apple Accelerate prefers threads=1. Pass threads=k to override; or set it yourself via BLAS.set_num_threads and call configure_blas!(exec; threads=BLAS.get_num_threads()) explicitly.

Returns the BLAS thread count actually set.

kspace_foreach!(body!, ks::AbstractMatrix, exec::Executor;
                scratch_factory = () -> nothing,
                progress = nothing)

Run body!(scratch, j, @view ks[:, j]) for j ∈ 1:size(ks, 2) under the given executor. scratch is task-local: scratch_factory() is called once per chunk (per thread / per worker), never per k. progress, if given, is a ProgressMeter.Progress — updates are funneled through a Channel so only one task touches the bar.

Returns nothing. Side effects on whatever body! mutates (bands/obs matrices, accumulated densities, etc.) are the caller's responsibility.

kspace_reduce!(body!, output, ks, exec; scratch_factory, scalar_init=nothing, progress=nothing)

Run body!(local_out, scratch, j, k) once per k-index. local_out is a per-task accumulator (one per chunk on threaded/distributed; just output itself on serial). After the loop, local_outs are summed into output via output .+= partial.

This is the right pattern for DOS, density, density-matrix, optical conductivity — anywhere each k contributes additively to a global accumulator. It is lock-free in the hot loop: each task writes only to its own local_out, and the master merges them once at the end.

progress, if given, is a ProgressMeter.Progress — updates are funneled through a Channel so only one task touches the bar.

output must support zero(output) and broadcasted .+=.

If scalar_init === nothing (default), kspace_reduce! returns output and the body's return value is ignored. If scalar_init is given (e.g. 0.0), each body return value is summed into a per-task accumulator and the total is returned in place of output — this is how to get a scalar side-channel (e.g. kinetic energy) under any executor, including distributed (where Ref capture in the body would be silently lost).

to_executor(x) -> Executor

Coerce a user-facing value (Executor, Symbol, or nothing) to an Executor. Symbols supported: :serial, :multithreaded/:threaded, :distributed, :auto. Falls back to SerialExec with a warning for unknown symbols (matches the previous multimode semantics).

padvec(v::AbstractVector, d::Int)

Make sure Vector v has length d, pad with zeros if needed.

spinorrotation(θ, n=[0,0,1])

Returns SU(2) spinor rotation matrix U parametrized by rotation angle θ and rotation axis n=[n1,n2,n3]. U = exp(-i (θ/2) σ.n).

Note: This routine will always normalize n.

scalar2vector(f0, N=3)

This is macro is a wrapper that takes as input a function f0(x::Float64) and adds a new dispatch f0(r1::Vector, r2::Vector) = f0(norm(r1-r2)) while making sure that r1 and r2 do not exceed length N.

getbands(H, ks::Union{DiscretePath, AbstractMatrix} [, As]; kwargs...)

Calculates the bands for operator H along discrete path ks and if operators As=[A1, A2, ...] are given, their expectaction values are calculated and stored for each eigenvector.

Note that ks is a discrete path object as returned by kpath(lat::Lattice,...).

Accepts the same keywords as geteigvals, geteigvecs, geteigen. In particular: format (:sparse or :dense) and num_bands::Int.

Returns a BandData object (with fields bands, obs, path).

Example

using LatticeQM

lat = Geometries.honeycomb()
h = Operators.graphene(lat)
ks = kpath(lat; num_points=200)
valley = Operators.valleyoperator(lat)

bands = getbands(h, ks, valley)

using Plots
plot(bands)

getdos(h, emin::Float64, emax::Float64, num=500; kwargs...)

Computes the density of states of operator h(k) on the entire Brillouine zone, discretized on a grid with $ k{lin} \times k{lin} $ points. and for the frequencies ωs=(ωmin, ω2, ..., ωmax). The paremter $\Gamma$ is the energy broadening.

Accepts the same kwargs as getdos(h, ωs; klin, Γ, kwargs...).

Note: the current implementation only works for a two-dimensional Brillouine zone. Might change in the future, but for now use dos(h, ks, ω; Γ) syntax if needed.

getdos(h, ωs; klin, Γ, kwargs...)

Computes the density of states of operator h(k) on the entire Brillouine zone, discretized on a grid with $ k{lin} \times k{lin} $ points. and for the frequencies ωs=(ω1, ω2, ...). The paremter $\Gamma$ is the energy broadening.

Accepts the same kwargs as dos(h, ks, ω).

Note: the current implementation only works for a two-dimensional Brillouine zone. Might change in the future, but for now use dos(h, ks, ω; Γ) syntax if needed.

getdos(h, ks, ω; Γ, parallel=true, format=:auto)

Computes the density of states of operator h(k) using the points ks=(k1,k2,...) and for the frequencies ω=(ω1, ω2, ...). The paremter $\Gamma$ is the energy broadening.

Mode can be :distributed or :serial, format can be :auto, :sparse or :dense.

opticalconductivity(frequencies, i, j, H, lat; kwargs...)

Compute the Kubo current-current correlator χ_ij(ω) for Hamiltonian H(k) using the Cartesian components i, j of the current operator defined on lat. Internally obtains J = getcurrentoperators(lat, H) and forwards to the low‑level routine, whose docstring gives the formula and the conversion from the returned χ(ω) to the optical conductivity σ(ω).

Common keywords: μ (chemical potential), Γ (broadening), T (temperature), plus diagonalization options forwarded to the eigen solver.

opticalconductivity(frequencies, H, J1, J2; klin, kwargs...)

Compute the Kubo current-current correlator χ(ω) on a regular klin × klin k‑grid using currents J1, J2. Convenience front‑end to the explicit‑grid method, whose docstring gives the formula and the χ(ω) → σ(ω) conversion.

opticalconductivity(frequencies, H, J1, J2, ks; μ=0.0, Γ=0.025, T=0.1, ...)

Evaluate the Kubo current-current correlator at the set of frequencies and k‑points ks (columns). H, J1, and J2 are callable with a k‑vector and return the Hamiltonian and current operators, respectively.

Returns a complex vector χ(ω) normalized by the number of k‑points N:

χ(ω) = -(i/N) Σ_k Σ_{m≠n} (f_m - f_n) ⟨m|J1|n⟩⟨n|J2†|m⟩ / (ϵ_m - ϵ_n - ω - iΓ)

where f = fermidirac(ϵ; μ, T) and |m⟩, ϵ_m are the Bloch eigenstates at k. This is not yet the optical conductivity. Obtain σ(ω) by subtracting the static (diamagnetic) piece and dividing by the photon energy:

σ(ω) = -(χ(ω) - χ(0)) / (ω + iΓ)

as done in test/test_linearresponse.jl and extra/examples/graphene/opticalconductivity.jl.

μ is the chemical potential, Γ the phenomenological broadening and T the temperature (all in the same energy units as H). Additional keyword arguments are forwarded to the eigenvalue solver.

solvehartreefock(h, v, ρ_init, filling; kwargs...)

Convenience wrapper around solveselfconsistent that constructs a HartreeFock functional from base Hamiltonian h and interaction kernel v. Returns the converged mean‑field solution and metadata.

solveselfconsistent(ρ0, mf::MeanfieldGenerator, filling, ks; kwargs...)
solveselfconsistent(ρ0, mf::MeanfieldGenerator, filling; klin, kwargs...)

Non‑mutating convenience wrappers around solveselfconsistent! that copy the initial density matrix ρ0, iterate the mean‑field functional mf (e.g., HartreeFock), and return the converged result together with energy and state.

The filling sets the target electronic filling (0–1 per spin). Supply either an explicit k‑grid ks or a grid resolution via klin (uses klin×klin).

Common keywords: iterations, tol, T (temperature), β (mixing), multimode (parallel), log_callback, residual_norm.

Returned energy decomposition

The returned mf carries the converged decomposition (physical-sign convention):

• mf.ϵkin — kinetic energy Tr[h₀ · ρ] (derived from the variational identity ϵkin = ϵband - 2(ϵH + ϵF), exact at the SCF fixed point)
• mf.ϵH — Hartree energy +½ nᵀ V₀ n (positive for repulsive V₀)
• mf.ϵF — Fock energy -½ Σ_L Σ_{ij} v_L[i,j] |ρ_L[i,j]|² (negative for repulsive v)
• mf.ϵband — band energy Tr[H_MF · ρ] = sum of occupied eigenvalues

The total HF energy (returned as the second tuple element) can be expressed in three equivalent ways at the SCF fixed point:

• E_HF = ϵkin + ϵH + ϵF (additive physical decomposition)
• E_HF = ϵband - ϵH - ϵF (band − double-counting)
• E_HF = ½ (ϵkin + ϵband) (½ Tr[ρ (h₀ + H_MF)])

When log_callback is supplied, the driver invokes it as log_callback(iter, ϵ, residual, info) if the callback can take a 4th arg (detected via applicable), with info = (ϵkin, ϵband, ϵH, ϵF, μ). Otherwise it falls back to the legacy 3-arg form log_callback(iter, ϵ, residual).

setfilling!(H, filling; nk=100, kwargs...)

Shift the chemical potential of Hamiltonian-like object H so that the resulting electronic filling matches filling (0–1 per spin). Internally uses Spectrum.chemicalpotential on a regular nk × nk grid and then applies addchemicalpotential!(H, -μ).

Keywords are forwarded to the underlying solver (e.g., temperature T, diagonalization format).

setfilling!(H, kgrid::Structure.Mesh, filling; kwargs...)

Variant that accepts an explicit k-grid kgrid (e.g., from regulargrid). Computes μ on that grid and shifts H in place.