API Reference

QuantumNumber{Names, Types}

A wrapper around NamedTuple for representing quantum numbers with cleaner display.

Type Parameters

• Names: Tuple of field names (Symbols)
• Types: Tuple type for values

Examples

# Create from keyword arguments
qn = QuantumNumber(site=1, spin=2)
qn.site  # 1
qn.spin  # 2

# Display: QN(site=1, spin=2)

Dof

Specification for a single degree of freedom (DOF).

Fields

• name::Symbol: Name of the DOF (e.g., :site, :spin, :orbital)
• size::Int: Number of states for this DOF
• labels::Union{Nothing, Vector}: Optional labels for states (e.g., [:up, :down])

SystemDofs{Q}

Defines all degrees of freedom in the system, with optional constraints, custom ordering, and symmetry blocks.

Type Parameter

• Q: QuantumNumber type with DOF names as keys

Fields

• dofs::Vector{Dof}: List of degrees of freedom
• valid_states::Vector{Q}: All valid quantum number combinations (sorted)
• blocks::Vector{UnitRange{Int}}: Index ranges for symmetry blocks (single block [1:N] if no blocking)
• qn_to_idx::Dict{Q, Int}: Fast O(1) lookup from quantum numbers to linear indices

Constructor

SystemDofs(dofs::Vector{Dof}; constraint=nothing, sortrule=collect(length(dofs):-1:1))

Arguments

• dofs: Vector of Dof specifications
• constraint: Optional function qn -> Bool to filter valid states (not saved after construction)
• sortrule: Vector or nested vector specifying sort priority and blocking (not saved after construction)
  • Default: [N, N-1, ..., 1] (reverse order: first dof varies fastest, no blocking)
  • Elements are indices into dofs array
  • Nested syntax for symmetry blocks:
    • [4, 3, 2, 1] - No blocking
    • [[4], 3, 2, 1] - Block by 4th DOF (e.g., spin)
    • [[4, 3], 2, 1] - Block by 4th DOF, then by 3rd within each block
    • [[4, 3, 2], 1] - Block by 4th, 3rd, 2nd sequentially

Design Rationale

• Default reverse order ensures first DOF (typically site) varies fastest for cache-friendly access
• Nested sortrule syntax naturally creates block-diagonal structure for conserved quantum numbers
• Blocks enable efficient matrix operations (multiplication, diagonalization) by exploiting sparsity

Examples

# Example 1: Default ordering (reverse, no explicit blocks)
dofs = SystemDofs([
    Dof(:site, 3),
    Dof(:spin, 2)
])
# sortrule = [2, 1]: spin slow, site fast
# valid_states = [QN(site=1,spin=1), QN(site=2,spin=1), QN(site=3,spin=1),
#                 QN(site=1,spin=2), QN(site=2,spin=2), QN(site=3,spin=2)]
# blocks = [1:6]  (single block containing all states)

# Example 2: With symmetry blocks (spin conservation)
dofs = SystemDofs([
    Dof(:site, 3),
    Dof(:spin, 2)
], sortrule = [[2], 1])
# blocks = [1:3, 4:6]  (spin-up block, spin-down block)

# Example 3: Nested blocks (spin and valley conservation)
dofs = SystemDofs([
    Dof(:site, 4),
    Dof(:orbital, 2),
    Dof(:spin, 2),
    Dof(:valley, 2)
], sortrule = [[4, 3], 2, 1])
# blocks = [1:8, 9:16, 17:24, 25:32]  (K↑, K↓, K'↑, K'↓)

# Example 4: With constraint
dofs = SystemDofs([
    Dof(:site, 4),
    Dof(:spin, 2)
], constraint = qn -> qn.site <= 2, sortrule = [[2], 1])
# blocks = [1:2, 3:4]  (only first 2 sites, with spin blocks)

Lattice{D, Q, T}

Represents a lattice structure mapping position states to real-space coordinates.

Type Parameters

• D: Spatial dimension (e.g., 2 for 2D, 3 for 3D)
• Q: QuantumNumber type for position states
• T: Coordinate element type (subtype of Real)

Fields

• position_dofs::Vector{Dof}: Degrees of freedom that define positions
• position_states::Vector{Q}: All valid position states
• coordinates::Vector{SVector{D,T}}: Coordinate for each position state (length D)
• vectors::Union{Nothing, Vector{SVector{D,T}}}: Supercell vectors for PBC (optional)

Notes

• position_states[i] corresponds to coordinates[i]
• vectors is set automatically when using tiling constructor

Examples

# 2D honeycomb lattice with 16 unit cells and 2 sublattices
position_dofs = [Dof(:moire_cell, 16), Dof(:sublattice, 2, [:A, :B])]
position_states = [QN(moire_cell=1, sublattice=1), QN(moire_cell=1, sublattice=2), ...]
coordinates = [[0.0, 0.0], [1.0, 0.0], ...]

lattice = Lattice(position_dofs, position_states, coordinates)

Bond{Q, T}

Represents a bond connecting one or more sites.

Fields

• states::Vector{Q}: Position states (length determines bond order: 1=onsite, 2=two-body, etc.)
• coordinates::Vector{SVector{D,T}}: Physical (unwrapped) coordinates for each site. For periodic bonds these may lie outside the simulation cell.
• icoordinates::Vector{SVector{D,T}}: Unit-cell origin for each site. icoordinates[n] is the lattice vector of the unit cell that contains coordinates[n]. For intra-cell bonds all entries are zero vectors. For bonds that cross a periodic boundary, the entry for the image site carries the lattice shift, e.g. [-2.0, 0.0].

Examples

# Onsite bond — icoordinates default to zero vectors
bond = Bond([state1], [coord1])

# Two-body intra-cell bond
bond = Bond([state1, state2], [coord1, coord2])

# Two-body bond crossing a periodic boundary (site 2 lives in cell [-2, 0])
bond = Bond([state1, state2], [coord1, coord2], [[0.0, 0.0], [-2.0, 0.0]])

FermionOp{Q} <: Operator

Represents a single fermion operator (creation or annihilation).

Fields

• qn::Q: Quantum number identifying the state
• dag::Bool: true = creation operator (c†), false = annihilation operator (c)

Operator

Abstract type for quantum operators. Subtypes include FermionOp (and potentially SpinOp, BosonOp, etc. in the future).

Operators{T}

Represents a product of operators with a coefficient.

Operators are stored in their original order (no automatic reordering). Reordering to canonical form (creation-annihilation alternating order: c†c c†c...) happens when building matrices/tensors via build_onebody_matrix or build_interaction_tensor.

Type Parameters

• T: Numeric type for coefficient (e.g., Float64, ComplexF64)

Fields

• value::T: Coefficient
• ops::Vector{<:Operator}: Operators in original order

Examples

i = QN(site=1, spin=1)
j = QN(site=2, spin=1)

# One-body: -t * c†_i c_j
op = Operators(-t, [cdag(i), c(j)])

# Two-body: U * c†_i↑ c_i↑ c†_i↓ c_i↓
op = Operators(U, [cdag(i_up), c(i_up), cdag(i_dn), c(i_dn)])

# Pair hopping: stored as-is, reordered when building tensor
op = Operators(J, [cdag(i_up), cdag(i_dn), c(j_dn), c(j_up)])

qn2linear(dofs::SystemDofs, qn) -> Int

Convert quantum numbers to linear index (position in valid_states). Uses O(1) hash table lookup for fast performance. Accepts QuantumNumber or NamedTuple.

Examples

dofs = SystemDofs([Dof(:site, 4), Dof(:spin, 2)])
idx = qn2linear(dofs, QN(site=1, spin=2))  # O(1) lookup
idx = qn2linear(dofs, (site=1, spin=2))    # also works

linear2qn(dofs::SystemDofs, idx::Int) -> QuantumNumber

Convert linear index back to quantum numbers.

Examples

dofs = site_spin_system(4)
qn = linear2qn(dofs, 2)  # returns QN(site=1, spin=2)

ndofs(dofs::SystemDofs) -> Int

Return number of degrees of freedom.

total_dim(dofs::SystemDofs) -> Int

Return total dimension (number of valid states).

site_spin_system(Nsite::Int) -> SystemDofs

Create a simple site+spin system (standard case).

get_coordinate(lattice::Lattice, state) -> Vector{T}

Return the coordinate for a given quantum state. The input can be QuantumNumber or NamedTuple, and can contain extra keys beyond the position DOFs.

Examples

lattice = Lattice(...)  # with position_dofs = [:moire_cell, :sublattice]

# Exact position state
coord = get_coordinate(lattice, QN(moire_cell=1, sublattice=2))

# Full quantum state with extra DOFs
coord = get_coordinate(lattice, QN(moire_cell=1, sublattice=2, spin=1, valley=1))

bonds(lattice, boundary, neighbors)

Generate bonds based on neighbor specification.

Arguments

• lattice::Lattice: The lattice structure (must have vectors set)
• boundary::NTuple{D, Symbol}: Boundary conditions, :p (periodic) or :o (open)
• neighbors: Neighbor specification:
  • Int: n-th nearest neighbor (0=onsite, 1=nearest, 2=next-nearest, ...)
  • Vector{Int}: multiple neighbor orders [n1, n2, ...]
  • Vector{<:Real}: specific distances [r1, r2, ...]

Returns

Vector{Bond} with appropriate number of sites per bond

Examples

# Create lattice with tiling (vectors auto-set)
lattice = Lattice(unitcell, (4, 4))

# Onsite bonds - returns Vector{Bond} with 1 site per bond
onsite = bonds(lattice, (:p, :p), 0)
# onsite[1] displays as: Bond{1}(QN(cell=1, sub=1) @ [0.000, 0.000])

# Nearest neighbor bonds - returns Vector{Bond} with 2 sites per bond
nn = bonds(lattice, (:p, :o), 1)
# nn[1] displays as: Bond{2}(QN(cell=1, sub=1) @ [0.000, 0.000], QN(cell=1, sub=2) @ [0.500, 0.289])

# Multiple neighbor orders
nn_and_nnn = bonds(lattice, (:p, :p), [1, 2])

# Specific distances
specific = bonds(lattice, (:p, :p), [1.0, 1.732])

is_positive_direction(delta::Vector) -> Bool

Check if delta vector is in "positive" direction using lexicographic ordering. The first non-zero component must be positive.

c(qn) -> FermionOp

Create an annihilation operator for quantum state qn.

cdag(qn) -> FermionOp

Create a creation operator (c†) for quantum state qn.

generate_onebody(dofs, bonds, value; order=(cdag, :i, c, :j), hc=true)

Generate one-body terms from bonds.

Arguments

• dofs::SystemDofs: Full DOF specification (position + internal DOFs)
• bonds::Vector{Bond}: One- or two-site bonds from bonds() function
• value: Coefficient, either:
  • Number: constant value for all internal DOF combinations
  • Function(delta, qn1, qn2) -> Number: custom function to constrain DOFs. delta = coord(op1_site) - coord(op2_site) in the user-specified order.

Keyword Arguments

• order::Tuple: Format (op_type1, label1, op_type2, label2) where labels are Symbols acting as site indices.
  • Default: (cdag, :i, c, :j) — standard hopping c†i cj
  • Same label (e.g. (cdag, :i, c, :i)) — on-site/chemical-potential term
  • Any operator order is accepted; operators are reordered to c†c alternating form automatically with the fermionic sign absorbed into the coefficient.
• hc::Bool: whether to include the Hermitian conjugate terms. Ignored (no hc added) when both labels are identical (on-site terms), since c†i ci is already self-contained and adding h.c. would double-count.

Returns

NamedTuple with three parallel Vector fields:

• .ops::Vector{Operators}: operator terms in c†c order, fermionic sign absorbed
• .delta::Vector{SVector{D,T}}: physical displacement coord(op1) - coord(op2) after reordering; sign is flipped for H.c. terms
• .irvec::Vector{SVector{D,T}}: unit-cell displacement after reordering; zero for intra-cell bonds, non-zero for bonds crossing a periodic boundary

NOTICE

• All internal DOFs are mixed unless constrained via the value function.

Examples

result = generate_onebody(dofs, nn_bonds, -1.0)
result.ops    # Vector{Operators}
result.delta  # physical displacements
result.irvec  # unit-cell displacements (use this for momentum-space HF)

# On-site chemical potential (no h.c. added automatically)
result = generate_onebody(dofs, onsite_bonds, -μ, order=(cdag, :i, c, :i))

# Unconventional input order — reordered internally with correct sign
result = generate_onebody(dofs, nn_bonds, 1.0, order=(c, :i, cdag, :j))

generate_twobody(dofs, bonds, value; order=(cdag, :i, c, :i, cdag, :j, c, :j))

Generate two-body interaction terms from bonds.

Arguments

• dofs::SystemDofs: Full DOF specification.
• bonds::Vector{Bond}: Bonds with 1–4 sites.
• value: Coefficient, either:
  • Number: constant for all internal DOF combinations.
  • Function(deltas, qn1, qn2, qn3, qn4) -> Number: custom amplitude function. deltas[m] = coord(op_m) - coord(op_4) for m = 1,2,3 (always 3 entries).

Keyword Arguments

• order::Tuple: Format (type1, label1, type2, label2, type3, label3, type4, label4) where labels are Symbols acting as free Einstein summation indices. Each unique Symbol is an independent position index summed over all injective assignments to bond sites (different Symbols → different sites).

  Common patterns:

  • (cdag, :i, c, :i, cdag, :j, c, :j) — density-density Σ{i≠j} ni n_j (default)
  • (cdag, :i, c, :i, cdag, :i, c, :i) — on-site Σi ni² (use with 1-site bonds)
  • (cdag, :i, c, :j, cdag, :i, c, :j) — exchange-type
  • (cdag, :i, c, :j, cdag, :i, c, :k) — three distinct positions (requires 3-site bond)

Returns

NamedTuple with three parallel Vector fields:

• .ops::Vector{Operators}: operator terms in creation-annihilation alternating order (c†c c†c), fermionic sign absorbed into the coefficient.
• .delta::Vector{NTuple{3, SVector{D,T}}}: physical displacements (δ1, δ2, δ3), δm = coord(op_m) - coord(op_4) after reordering.
• .irvec::Vector{NTuple{3, SVector{D,T}}}: unit-cell displacements (τ1, τ2, τ3), τm = icoord(op_m) - icoord(op_4) after reordering.

Examples

# Nearest-neighbor Coulomb: V Σ_{i≠j,σσ'} n_{iσ} n_{jσ'}  (full ordered-pair sum)
twobody = generate_twobody(dofs, nn_bonds, V)

# Onsite Hubbard U: Σ_i U n_{i↑} n_{i↓}
twobody = generate_twobody(dofs, onsite_bonds,
    (deltas, qn1, qn2, qn3, qn4) ->
        (qn1.spin == qn2.spin) && (qn3.spin == qn4.spin) && (qn1.spin !== qn3.spin) ? U/2 : 0.0,
    order = (cdag, :i, c, :i, cdag, :i, c, :i))

# NN Coulomb in x-direction only (filter by deltas[1] = coord(i) - coord(j))
twobody = generate_twobody(dofs, nn_bonds,
    (deltas, qn1, qn2, qn3, qn4) -> abs(deltas[1][1]) ≈ 1.0 && iszero(deltas[1][2]) ? V : 0.0)

# Pair hopping: Σ_{i≠j} J c†_{i↑} c†_{i↓} c_{j↓} c_{j↑}
twobody = generate_twobody(dofs, nn_bonds,
    (deltas, qn1, qn2, qn3, qn4) ->
        (qn1.spin, qn2.spin, qn3.spin, qn4.spin) == (1, 2, 2, 1) ? J : 0.0,
    order = (cdag, :i, cdag, :i, c, :j, c, :j))

build_onebody_matrix(dofs::SystemDofs, ops::AbstractVector{<:Operators})

Build Hermitian one-body Hamiltonian matrix from Operators. Each Operators must have exactly 2 FermionOp. Operators are reordered to c†c format internally.

build_interaction_tensor(dofs::SystemDofs, ops::AbstractVector{<:Operators})

Build two-body interaction tensor from Operators. Each Operators must have exactly 4 FermionOp. Operators are reordered to creation-annihilation alternating order (c†c c†c) internally.

The tensor V[i,j,k,l] corresponds to c†i cj c†k cl in creation-annihilation alternating order.

build_T(dofs, ops)

Build sparse one-body Hamiltonian (N×N) from one-body operators. Only stores elements within the same symmetry block (from dofs.blocks). Returns a sparse Hermitian matrix.

build_U(dofs, ops; include_fock=true)

Build sparse HF interaction matrix (N²×N²) from two-body operators. Applies the antisymmetrization formula directly, skipping the V tensor. Only stores columns where (k,l) share a symmetry block. Returns a sparse matrix of size (N²×N²).

Set include_fock=false for Hartree-only (no exchange) calculation. The 4-term decomposition is:

• Hartree: U_{ij,kl} += V_{ijkl} + V_{klij}
• Fock: U_{ij,kl} -= V_{kjil} + V_{ilkj} (skipped if include_fock=false)

solve_hfr(dofs, ops, block_occupations; kwargs...)

Solve Hartree-Fock equations using self-consistent field (SCF) iteration.

When the interaction is strong (large U/t), the HF energy landscape has multiple local minima. Use n_restarts > 1 to run multiple random initializations and automatically select the lowest-energy converged solution.

Arguments

• dofs::SystemDofs: System degrees of freedom
• ops: All operators (one-body with 2 FermionOp and two-body with 4 FermionOp)
• block_occupations: Particle number — Int for single block, Vector{Int} for multiple blocks

Keyword Arguments

• temperature::Float64 = 0.0: Temperature (0 for ground state)
• max_iter::Int = 1000: Maximum SCF iterations per restart
• tol::Float64 = 1e-6: Convergence tolerance for Green's function residual
• mix_alpha::Float64 = 0.5: Linear mixing parameter (0 < α ≤ 1), used when DIIS history is insufficient or disabled.
• diis_m::Int = 8: DIIS history window length. Stores the last diis_m iterates and solves a small linear system to extrapolate the optimal density matrix. Set to 0 to disable DIIS and use pure linear mixing.
• G_init = nothing: Initial Green's function for the first restart
• ene_cutoff::Float64 = 100.0: Energy cutoff for exp overflow at finite T
• n_restarts::Int = 1: Number of random restarts; returns lowest-energy converged result
• seed::Union{Nothing, Int} = nothing: Random seed for reproducibility
• include_fock::Bool = true: Include Fock exchange terms. Set to false for Hartree-only.
• verbose::Bool = true: Print iteration information
• field_strength::Float64 = 0.0: Amplitude of a random symmetry-breaking field added to h_eff for the first n_warmup iterations of each restart (linearly decayed to zero). Set to a value comparable to the hopping scale to escape paramagnetic or other high-symmetry saddle points without prior knowledge of the ordered phase.
• n_warmup::Int = 15: Number of warmup iterations during which the symmetry-breaking field is active. Should be well below the typical SCF convergence count.

Returns

NamedTuple: (G, eigenvalues, eigenvectors, energies, mu_list, converged, iterations, residual, ncond, sz)

• ncond: Total particle number Σ G[i,i]
• sz: Total spin Sz (only if dofs contains a :spin Dof with size=2, otherwise nothing)

build_Tr(dofs, ops, irvec) -> NamedTuple

Parse one-body operators and build the real-space hopping table T_{ab}(r).

Scans irvec to find all unique displacement vectors, then accumulates hopping amplitudes into one d_int × d_int matrix per displacement. Non-one-body operators (i.e. those with length(op.ops) ≠ 2) in ops are silently skipped, so the full operator list from generate_onebody (which may include metadata entries) can be passed directly.

Arguments

• dofs::SystemDofs: DOFs of one magnetic unit cell.
• ops::Vector{<:Operators}: operator terms; only those with exactly 2 FermionOps (one-body) are used. ops may contain operators of other sizes — they are ignored.
• irvec::Vector{<:AbstractVector{<:Real}}: displacement vector per operator term, paired one-to-one with ops.

Returns

NamedTuple (mats, rvec) where:

• mats::Vector{Matrix}: one d_int × d_int hopping matrix per unique displacement; element type matches eltype of the operator values (e.g. Float64 or ComplexF64).
• rvec::Vector{Vector{Float64}}: the corresponding displacement vectors, same order as mats.

If ops is empty, emits a warning and returns (mats=Matrix{Float64}[], rvec=Vector{Float64}[]).

build_Tk(T_r) -> Function

Return a closure that evaluates the kinetic Hamiltonian at a given k-point:

T_{ab}(k) = Σ_r  T_{ab}(r) · exp(i k · r)

Arguments

• T_r: Real-space hopping table from build_Tr (NamedTuple with mats and delta).

Returns

A callable T_func(k) returning Matrix{ComplexF64} of shape (d_int, d_int), where k is a momentum vector. Returns nothing if T_r.mats is empty (no hopping terms); the caller should skip the kinetic contribution in that case.

build_Vr(dofs, ops, irvec) -> NamedTuple

Parse two-body operators and build the real-space interaction table V̄^{abcd}(τ1, τ2, τ3).

Scans irvec for unique (τ1,τ2,τ3) triples, then accumulates interaction amplitudes into one d_int×d_int×d_int×d_int array per unique triple. Non-four-body operators in ops are silently skipped.

Arguments

• dofs::SystemDofs: DOFs of one magnetic unit cell; d_int = length(dofs.valid_states).
• ops::Vector{<:Operators}: Two-body operators in creation-annihilation alternating order (4 FermionOp entries); typically generate_twobody(...).ops.
• irvec: Per-operator unit-cell displacements (τ1, τ2, τ3) paired one-to-one with ops; typically generate_twobody(...).irvec.

Returns

NamedTuple (mats, taus) where:

• mats::Vector{Array{T,4}}: one d_int×d_int×d_int×d_int array per unique (τ1,τ2,τ3); mats[n][a,b,c,d] = V̄^{abcd} for the n-th displacement triple.
• taus::Vector{NTuple{3,Vector{Float64}}}: the corresponding (τ1,τ2,τ3) triples, same order as mats.

Notes

• Raw V̄ values are stored without antisymmetrization; the HF self-energy symmetrization is handled inside build_heff_k!.
• The density-density special case (all τ1=τ3=0) is detected downstream in build_heff_k! to select the real-space path.

build_Vk(V_r, kgrid) -> Function

Fourier-transform V̄^{abcd}(τ1, τ2, τ3) to the three-momentum interaction kernel:

Ṽ^{abcd}(k1, k2, k3) = Σ_{τ1,τ2,τ3}  V̄^{abcd}(τ1, τ2, τ3)
                        · exp(-i k1·τ1 + i k2·τ2 - i k3·τ3)

Returns a closure (k1_idx, k2_idx, k3_idx) -> Array{ComplexF64, 4} that evaluates the interaction kernel at any three k-points on demand. The fourth momentum k4 = k1 + k3 - k2 is fixed by momentum conservation.

Arguments

• V_r: Real-space interaction table from build_Vr.
• kgrid: k-grid struct providing k_points and nk.

Returns

A callable Ṽ(k1_idx, k2_idx, k3_idx) returning Array{ComplexF64, 4} of shape (d_int, d_int, d_int, d_int).

Notes

• Only needed for the general (non-density-density) path in build_heff_k!.
• For density-density interactions (all τ1=τ2, τ3=0 in Vr.entries), `buildheffk!uses real-space directly fromVr` without calling this.

build_Uk(V_k) -> Function

Return a closure that, given two k-points k and q, assembles the antisymmetrized HF interaction matrix of shape (d²,d²):

U[α+(β-1)d, μ+(ν-1)d] =  Ṽ^{μναβ}(k,k,q)   [Hartree 1]
                         + Ṽ^{αβμν}(q,q,k)   [Hartree 2]
                         - Ṽ^{μβαν}(k,q,q)   [Fock 1]
                         - Ṽ^{ανμβ}(q,k,k)   [Fock 2]

This is the k-space analogue of the real-space build_U, encoding all four Wick contraction channels (Theory §4). The row index α+(β-1)d corresponds to Σ^{αβ}(q) and the column index μ+(ν-1)d corresponds to G^{μν}(k), both following Julia's column-major convention (i.e. vec(M)[α+(β-1)*d] == M[α,β]).

Once U_func = build_Uk(V_k), the HF self-energy is

Σ(q) = (1/Nk) Σ_k reshape(U_func(k, q) * vec(G(k)), d, d)

Returns nothing if V_k is nothing.

build_kpoints(unitcell_vectors, box_size) -> Vector{Vector{Float64}}

Generate the k-point grid for a D-dimensional lattice.

k = (m1/n1)*b1 + (m2/n2)*b2 + ...  for m_i in 0:n_i-1

where bi are reciprocal lattice vectors satisfying ai · bj = 2π δ{ij}.

Arguments

• unitcell_vectors::Vector{Vector{Float64}}: Direct lattice vectors [a1, a2, ...].
• box_size::NTuple{D,Int}: Number of unit cells along each direction.

initialize_green_k(Nk, d; G_init=nothing, n_fill=nothing, rng=default_rng()) -> Array{ComplexF64, 3}

Initialize the k-space Green's function stored as G[a, b, k_idx] of shape (d, d, Nk) (column-major: each G[:,:,ki] is a contiguous d×d matrix in memory).

If G_init is provided, validates shape and Hermiticity at each k and returns it.

Otherwise, each G(k) is a random Hermitian matrix with eigenvalues in [0, 1] and mean filling n_fill / (Nk * d) per orbital. This initialization has physical magnitude (unlike a small perturbation around zero), so the first SCF step produces a non-trivial Hartree–Fock field that can break symmetry and escape paramagnetic saddle points. If n_fill is not provided, eigenvalues are drawn uniformly from [0, 1].

Arguments

• Nk::Int: Number of k-points
• d::Int: Internal dimension per unit cell

Keyword Arguments

• G_init: Pre-initialized G array of shape (d, d, Nk), or nothing
• n_fill::Union{Nothing,Int}: Total electron number; used to set the mean occupation per k-point so that the initial filling is approximately correct.
• rng: Random number generator

Returns

Array{ComplexF64, 3} of shape (d, d, Nk).

green_k_to_tau(G_k, kpoints, taus) -> Vector{Matrix{ComplexF64}}

Compute G(τ) at each displacement in taus via direct Fourier sum:

G^{ab}(τ) = (1/Nk) Σ_k  G^{ab}(k) · exp(-i k · τ)

Arguments

• G_k::Array{ComplexF64, 3}: Shape (d, d, Nk).
• kpoints::AbstractVector{<:AbstractVector{Float64}}: k-point list of length Nk.
• taus::Vector{Vector{Float64}}: displacement vectors at which to evaluate G(r).

solve_hfk(dofs, onebody, twobody, kpoints, n_electrons; kwargs...)

Solve Hartree-Fock in momentum space via SCF iteration.

Each k-point is diagonalized independently (parallelized over Threads.nthreads() threads). The magnetic unit cell is encoded in dofs; kpoints samples its Brillouin zone. For broken-symmetry phases (AFM, CDW, …) pass an enlarged magnetic unit cell.

Arguments

• dofs::SystemDofs: Internal DOFs of one magnetic unit cell.
• onebody: Result of generate_onebody(...); needs .ops and .irvec.
• twobody: Result of generate_twobody(...); needs .ops and .irvec.
• kpoints::AbstractVector{<:AbstractVector{Float64}}: k-points (use build_kpoints for uniform grids).
• n_electrons::Int: Total electron number across all k-points.

Keyword Arguments

• temperature::Float64 = 0.0
• max_iter::Int = 1000
• tol::Float64 = 1e-8: Convergence threshold ‖ΔG‖_F / (Nk·d²).
• diis_m::Int = 8: DIIS history window (0 = disabled).
• G_init = nothing: Initial G_k of shape (d, d, Nk); random if nothing.
• ene_cutoff::Float64 = 100.0
• n_restarts::Int = 1: Random restarts; returns lowest-energy converged result.
• seed::Union{Nothing,Int} = nothing
• include_fock::Bool = true
• verbose::Bool = true
• field_strength::Float64 = 0.0: Amplitude of the random symmetry-breaking field added to Hk for the first `nwarmupiterations of each restart. The field is a random Hermitian matrix (same for all k, independent per restart), linearly decayed to zero over the warmup phase. Set to a value comparable to the hopping scale (e.g.1.0) combined withn_restarts > 1` to escape paramagnetic saddle points without prior knowledge of the ordered phase.
• n_warmup::Int = 15: Number of warmup iterations during which the symmetry-breaking field is active (linearly decayed to zero). Only used when field_strength > 0. Should be well below the typical SCF convergence count so the system can converge freely after the field has vanished.

Returns

NamedTuple: G_k, eigenvalues, eigenvectors, energies, mu, kpoints, converged, iterations, residual, ncond.

energy_bands(dofs, onebody, twobody, kgrid, G_k, qpoints; include_fock=true)
    -> (bands::Matrix{Float64}, eigenvectors::Array{ComplexF64,3})

Compute strict HF band energies along arbitrary q-points from a converged G_k on a uniform k-grid. Returns the band matrix (d × Nq) and the eigenvector tensor (d × d × Nq), where eigenvectors[:, n, qi] is the n-th eigenvector at q-point qi (columns sorted by ascending eigenvalue).

local_spin(dofs, G; spin_dof=:spin, kwargs...) -> Vector{NamedTuple}

Compute the local magnetic moment on each site from a converged HF density matrix, with optional label-based filtering via keyword arguments.

A site is defined by all quantum numbers except the spin DOF, so the function works for any system (simple site+spin, multi-orbital, supercell, …).

Compatible with both HF solvers:

• HFr: pass result.G — a real-space density matrix of shape (N, N)
• HFk: pass result.G_k — a k-space density matrix of shape (d, d, Nk); the BZ average G_loc = Σ_k G_k / Nk is computed automatically.

Requirements

dofs must contain exactly one spin DOF of size 2 (default name :spin). Spin label convention: value 1 → spin-up, value 2 → spin-down.

Keyword arguments

• spin_dof: name of the spin DOF (default :spin)
• any other keyword: filter criterion on the site label. The value can be a single integer or a vector of integers. Only sites matching all criteria are returned.

Returns

A Vector of NamedTuples, one per (matching) site, ordered by each site's first linear index in dofs.valid_states. Each entry contains:

Spin-component formulas

Given G_loc[a,b] = ⟨c†_a c_b⟩:

• mz = (G[↑,↑] − G[↓,↓]) / 2
• mx = Re(G[↑,↓])
• my = Im(G[↑,↓])

Examples

result = solve_hfk(...)

# All sites
mags = local_spin(dofs, result.G_k)

# Only orbital=1 (a orbital) across all sites
mags_a = local_spin(dofs, result.G_k; orbital=1)

# Only B sublattice
mags_B = local_spin(dofs, result.G_k; sub=2)

# B sublattice at cells 3 and 4
mags_B34 = local_spin(dofs, result.G_k; sub=2, cell=[3,4])

print_spin(mags_B)

print_spin(mags; digits=4, io=stdout)

Pretty-print the output of local_spin.

Columns: site label | n | mx | my | mz | |m| | θ(°) | φ(°)

Example output

site                     n       mx       my       mz      |m|     θ(°)     φ(°)
────────────────────────────────────────────────────────────────────────────────
1                   1.0000   0.0000   0.0000   0.2341   0.2341     0.00     0.00
2                   1.0000   0.0000   0.0000  -0.2341   0.2341   180.00     0.00

spin_structure_factor(mags, lattice, qpoints) -> Vector{Float64}

Compute the magnetic structure factor

\[S(\mathbf{q}) = \frac{1}{N_s} \left|\sum_s e^{i\mathbf{q}\cdot\mathbf{d}_s}\, \vec{m}_s\right|^2\]

at each q-point, where the sum runs over all sites, \mathbf{d}_s is the real-space coordinate of site s (looked up from lattice), and \vec{m}_s = (m_x, m_y, m_z) is its local magnetic moment (e.g. from local_spin).

Arguments

• mags: output of local_spin
• lattice: the Lattice used to build the system; site coordinates are extracted automatically from mag.label via get_coordinate
• qpoints: collection of q-vectors to evaluate at (same format as build_kpoints)

Returns

Vector{Float64} of length length(qpoints).

Example

mags = local_spin(dofs, result.G_k)
qpts = build_kpoints([a1, a2], (100, 100))
Sq   = spin_structure_factor(mags, lattice, qpts)

magnetic_ordering_wavevector(mags, lattice, qpoints) -> NamedTuple

Find the magnetic ordering wavevector Q as the q-point that maximises spin_structure_factor.

Arguments

• mags: output of local_spin
• lattice: the Lattice used to build the system
• qpoints: collection of q-vectors to search over (same format as build_kpoints)

Returns

A NamedTuple with fields:

Example

mags = local_spin(dofs, result.G_k)
qpts = build_kpoints([a1, a2], (100, 100))
res  = magnetic_ordering_wavevector(mags, lattice, qpts)
println("Q = ", res.Q, "  S(Q) = ", res.Sq_max)

local_charge(dofs, G; kwargs...) -> Vector{NamedTuple}

Compute the local charge (occupation) for each valid quantum state from a converged HF density matrix.

All DOFs (spin, orbital, sublattice, …) are treated equally: no DOF is singled out or summed over. Each valid state in dofs gets its own occupation n = G_loc[i, i]. The full quantum-number label is preserved so that charge_structure_factor can extract spatial coordinates via the Lattice.

Compatible with both HF solvers:

• HFk: pass result.G_k — BZ average is computed automatically.
• HFr: pass result.G — used directly.

Keyword arguments

Any keyword filters on the quantum-number label (same convention as local_spin).

Returns

Vector of NamedTuples with fields label and n, one per matching state, sorted by linear index.

Example

charges  = local_charge(dofs, result.G_k)
charges_A = local_charge(dofs, result.G_k; sub=1)

charge_structure_factor(charges, lattice, qpoints) -> Vector{Float64}

Compute the charge structure factor

\[N(\mathbf{q}) = \frac{1}{N_s} \left|\sum_s e^{i\mathbf{q}\cdot\mathbf{d}_s}\, \delta n_s\right|^2\]

where $\delta n_s = n_s - \bar{n}$ is the deviation of the local charge from its spatial average. Subtracting the mean suppresses the trivial $\mathbf{q}=0$ peak from uniform filling, revealing true CDW modulations.

Arguments

• charges: output of local_charge
• lattice: the Lattice used to build the system; site coordinates are extracted automatically via get_coordinate
• qpoints: collection of q-vectors to evaluate at (same format as build_kpoints)

Returns

Vector{Float64} of length length(qpoints).

Example

charges = local_charge(dofs, result.G_k)
qpts    = build_kpoints([a1, a2], (100, 100))
Nq      = charge_structure_factor(charges, lattice, qpts)

charge_ordering_wavevector(charges, lattice, qpoints) -> NamedTuple

Find the charge ordering wavevector Q as the q-point that maximises charge_structure_factor.

Returns

A NamedTuple with fields:

Example

charges = local_charge(dofs, result.G_k)
qpts    = build_kpoints([a1, a2], (100, 100))
res     = charge_ordering_wavevector(charges, lattice, qpts)
println("Q = ", res.Q, "  N(Q) = ", res.Nq_max)

plot_lattice(lattice; kwargs...) -> Figure

Visualize a Lattice: sites (optionally colored by a DOF) and one or more sets of bonds drawn as line segments.

Requires Makie to be loaded before calling:

using GLMakie    # interactive
using CairoMakie # save to file

Arguments

• lattice: the Lattice to visualize

Keyword arguments

Example

using CairoMakie
unitcell = Lattice([Dof(:cell,1), Dof(:sub,2,[:A,:B])],
                   [QN(cell=1,sub=1), QN(cell=1,sub=2)],
                   [[0.0,0.0],[1.0,0.0]]; vectors=[[0.0,√3],[1.5,√3/2]])
lattice  = Lattice(unitcell, (4, 4))
nn  = bonds(lattice, (:p,:p), 1)
nnn = bonds(lattice, (:p,:p), 2)

fig = plot_lattice(lattice;
    bond_lists   = [nn, nnn],
    bond_colors  = [:tomato, :steelblue],
    bond_labels  = ["NN", "NNN"],
    site_groupby = :sub,
    site_colors  = [:seagreen, :mediumpurple])
save("lattice.png", fig)

plot_spin(mags, positions; kwargs...) -> Figure

Visualize local magnetic moments as arrows on a lattice.

Requires Makie to be loaded before calling:

using GLMakie    # interactive 3D
using CairoMakie # save to file
using WGLMakie   # Jupyter notebook

Arguments

• mags: output of local_spin
• positions: coordinates of each site, same order as mags. Each element is a 2- or 3-component vector [x, y] or [x, y, z].

Keyword arguments

Layout

Single top-view (xy) panel. All three spin components are encoded simultaneously:

• Arrow: direction = (mx, my), length ∝ sqrt(mx²+my²) (in-plane projection)
• ⊙ dot: mz > 0 (spin out of page), size ∝ mz
• ⊗ cross: mz < 0 (spin into page), size ∝ |mz|

Sites with purely in-plane spins show only arrows; purely z-polarized sites show only ⊙/⊗; canted spins show both simultaneously.

Example

using CairoMakie
result = solve_hfk(...)
mags   = local_spin(dofs, result.G_k)

nn = bonds(lattice, (:p, :p), 1)
fig = plot_spin(mags, unitcell.coordinates;
    bonds         = nn,
    unitcell_vecs = [[3/2, √3/2], [0.0, √3]],
    title         = "KMH ground state")
save("magnetization.pdf", fig)

solve_ph_excitations(dofs, onebody, twobody, hf_result, qpoints,
                     reciprocal_vecs; solver=:TDA, n_list=nothing,
                     eta=1e-10, verbose=true) -> NamedTuple

Compute the particle-hole excitation spectrum along a q-path.

Shifted momenta (k+q, p-q, …) are folded into the first BZ. If the folded point falls on the HF k-mesh, its eigensystem is reused; otherwise it is recomputed on the fly via energy_bands.

Particle-hole pairs are auto-detected per k-point from eigenvalues vs mu.

Keyword n_list

Optional vector of band indices allowed as particle (unoccupied) states. Default nothing means all bands are candidates. Use e.g. n_list=[1,2] to restrict the particle-hole space to low-energy bands only.

Keyword eta

Regularization added to the diagonal of the Bosonic BdG matrix before Cholesky decomposition (RPA only). Default 1e-10. Increase (e.g. 1e-6) for systems where the matrix is nearly singular near Goldstone modes.

Keyword solver

• :TDA (default): Tamm-Dancoff (Hermitian, M×M). Builds only A.
• :RPA: Full RPA. Builds A, B, and D, then solves the 2M×2M eigenvalue problem: (A B; -B† D)(X; Y) = ε(X; Y). The D block uses D(q) = -A(-q)* when the ±q occupation spaces match, otherwise it is computed directly with the residual occupation factor.

rd(x) -> Number

Round a number to PRECISION decimal places. Works with Real and Complex.

Global precision for numerical rounding (number of decimal places). All floating-point values (coordinates, hopping amplitudes, interactions, phases, etc.) are rounded to this precision to avoid floating-point comparison issues.