Skip to content

Wannierization Workflows

Wannierization is a controlled passage from validated Bloch source data to a localized representation, not a cosmetic replacement of a band plot by a set of orbitals. The calculation must declare which subspace is retained, whether that subspace was fixed by spectral isolation or selected from entangled bands, which gauge was chosen inside it, which operators were transformed, and what independent evidence licenses the resulting files.

This page owns the numerical selection, construction, convergence, validation, and reproducibility workflow. It stops before interpreting the localized Hamiltonian as a material model, classifying a topological phase, assigning a polarization, choosing screened interactions, or computing a transport coefficient. Those claims require their canonical physical owners.

Required background. Wannier Functions supplies projectors, frame gauge, the Bloch–Wannier transform, localization, centers, and obstruction theory. Band Structure Workflows supplies the fixed crystal and state, converged density, full-zone source states and overlaps, symmetry data, operator matrices, and immutable provenance that this workflow consumes.

Begin with a sentence that can fail:

For this fixed Bloch source, construct this rank-JJ localized representation, reproduce these energies and operators over this domain to these tolerances, and archive enough evidence to license this downstream use.

The sentence separates three objects that are often blurred together:

  1. the selected projector PkP_{\mathbf k};
  2. a localized frame spanning that projector; and
  3. the finite-range Hamiltonian and requested operator matrices exported from that frame.

For a spectrally isolated composite manifold, the first object is already fixed by the source Hamiltonian. Wannierization then chooses a convenient frame inside it. For entangled bands, selecting a rank-JJ projector from a larger source space is an additional model choice. A smooth, compact-looking orbital cannot erase that distinction.

The target accuracy must be attached to a downstream use. A representation acceptable for a density-of-states plot may fail a velocity-matrix tolerance needed for transport. A Hamiltonian that reproduces energies can still give wrong spin, dipole, or Berry-related matrix elements. Conversely, an overly strict tolerance on an object that will never be used wastes computation without strengthening the claim.

A complete workflow therefore does five things in order:

  1. freezes and audits the source physical problem;
  2. declares rank, windows, trials, symmetry, and operator requirements;
  3. selects a projector when necessary and optimizes a gauge inside it;
  4. validates energies, projectors, operators, and real-space truncation on withheld full-zone data; and
  5. archives the accepted and failed alternatives with an explicit stopping rule.

If admissible windows or trial orbitals yield materially different projectors or operators, the result is representation instability. More optimizer steps do not convert that instability into convergence.

Use the same ten labels as the Computational Quantum Matter gateway. Fill every field before a production export; write unresolved when a required choice or datum is not known.

  1. Physical problem. State the material or model, degrees of freedom, dimension, geometry, boundaries, specimen or unit cell, and preparation.
  2. State and limit. Give ensemble, filling or chemical potential, temperature, field, disorder, drive, equilibrium status, finite-size or thermodynamic target, and the required order of limits.
  3. Claim and accuracy. Name the target observable or decision, units, normalization, resolution, tolerance, and motivating physical comparison.
  4. Representation and provenance. Record the Hamiltonian or functional, active and outer source spaces, source metric and basis, spinor and orbital conventions, material parameters, artifact identifiers, omitted degrees of freedom, and parameter sources.
  5. Method and controlled domain. State the isolated-projector or disentangled-subspace branch, projection and localization method, symmetry constraints, controlled approximations, and why the construction can estimate the requested object.
  6. Finite numerical problem. Declare cell and boundary conventions, source and holdout meshes, rank, outer and frozen windows when applicable, trial orbitals, real-space cutoff, tolerances, starts or seeds, and stored source bands.
  7. Estimator and forward model. Define spread, projector, interpolation, tail, symmetry, and requested operator metrics and the map from the localized representation to the downstream quantity.
  8. Convergence and uncertainty. Separate source-model, rank, window, trial, mesh, real-space truncation, optimizer, symmetry, operator, fit, and downstream-forward-model uncertainties, including covariance where it matters.
  9. Verification, validation, and provenance. List orthonormality, Hermiticity, frozen containment, symmetry sewing, exact limits, independent full-zone holdout, alternative windows and trials, raw artifacts, versions, hashes, and failed runs.
  10. Licensed claim and stopping rule. State the strongest supported representation claim, credible alternatives, missing physics, falsifier, stop condition, and next canonical owner.

The record is prospective: tolerances, comparison metrics, and stopping rules must be fixed before inspecting the preferred result. A final band overlay is not a workflow record, and the smallest-spread run is not automatically the accepted one.

Wannierization begins only after the source-data owner has frozen the physical problem. A changed density, structure, functional, pseudopotential, magnetic state, relativistic convention, occupation rule, or stored band set is a new source problem rather than a Wannierization control.

Let the source outer frame Ψk\Psi_{\mathbf k} be an isometry from an Nout(k)N_{\rm out}(\mathbf k)-dimensional coefficient space to the physical Bloch fiber:

Ψk†Ψk=INout(k).\Psi_{\mathbf k}^{\dagger}\Psi_{\mathbf k} = I_{N_{\rm out}(\mathbf k)}.

In a nonorthogonal coefficient basis, archive the physical overlap matrix and use, at each momentum,

Ck†SkCk=I.C_{\mathbf k}^{\dagger} S_{\mathbf k} C_{\mathbf k} = I.

Same-momentum normalization and trial projections use SkS_{\mathbf k}. A neighbor link instead needs the archived cross-momentum physical overlap and sewing map. In coefficient notation its generic form is

M(k,b)=Ck†Sk,k+bCk+b.M(\mathbf k,\mathbf b) = C_{\mathbf k}^{\dagger} S_{\mathbf k,\mathbf k+\mathbf b} C_{\mathbf k+\mathbf b}.

A projected physical operator uses Ck†OkphysCkC_{\mathbf k}^{\dagger}O^{\rm phys}_{\mathbf k}C_{\mathbf k}. An expression such as C†SOcoeffCC^{\dagger}S O^{\rm coeff}C is valid only after declaring the coefficient-action convention that defines OcoeffO^{\rm coeff}. Euclidean normalization of generalized eigenvectors changes the physical subspace and invalidates later singular-value tests.

The immutable handoff should contain at least:

  • lattice vectors, origin, handedness, orbital embedding, reciprocal basis, and boundary-sewing convention;
  • electron count, filling or chemical potential, physical temperature, magnetic and spinor conventions, and the accepted density or represented Hamiltonian;
  • source eigenvalues, physical states or metric-aware coefficients, stored band count, mesh points and weights, and symmetry maps;
  • overlap matrices on every neighbor link, including reciprocal-boundary sewing phases;
  • physical source matrices for every requested operator, including nonlocal or augmentation contributions where applicable;
  • full-zone isolation or outer-space coverage evidence; and
  • code, environment, input, pseudopotential, raw-output, and artifact hashes.

An isolated-manifold request needs the minimum direct separation from the complement over the full Brillouin zone. An entangled request instead needs a source outer space wide enough to support every tested rank and frozen-state constraint. A high-symmetry path establishes neither condition.

The source handoff also fixes conventions that cannot be reconstructed from eigenvalues alone: whether spinors are explicit, how Kramers partners are ordered, how orbitals cross a reciprocal boundary, how the cell origin and site axes are embedded, and which velocity or current operator the source actually evaluated. If any of these are unresolved, stop before optimizing a gauge.

Separate Isolated Manifolds from Disentangled Subspaces

Section titled “Separate Isolated Manifolds from Disentangled Subspaces”

Choose a rank-JJ subspace with a semiunitary matrix VkV_{\mathbf k} and a frame inside it with UkU_{\mathbf k}:

Vk†Vk=IJ,Uk∈U(J),Qk=VkUk.V_{\mathbf k}^{\dagger}V_{\mathbf k}=I_J, \qquad U_{\mathbf k}\in U(J), \qquad Q_{\mathbf k}=V_{\mathbf k}U_{\mathbf k}.

The selected physical frame and projector are

Φk=ΨkQk,\Phi_{\mathbf k} = \Psi_{\mathbf k}Q_{\mathbf k}, Pk=ΨkVkVk†Ψk†.P_{\mathbf k} = \Psi_{\mathbf k} V_{\mathbf k}V_{\mathbf k}^{\dagger} \Psi_{\mathbf k}^{\dagger}.

A change Uk↦UkWkU_{\mathbf k}\mapsto U_{\mathbf k}W_{\mathbf k} with Wk∈U(J)W_{\mathbf k}\in U(J) leaves PkP_{\mathbf k} fixed. It changes orbital centers, spreads, onsite matrices, and hoppings, but it is an exact basis change when every state and operator is transformed consistently.

For a genuinely isolated manifold, the rank-JJ spectral projector is directly separated from excluded states at every momentum. Internal crossings and degeneracies are allowed. In this branch Nout(k)=JN_{\rm out}(\mathbf k)=J and there is no disentanglement, outer-window optimization, or frozen-window optimization. VkV_{\mathbf k} merely represents the already fixed projector.

Direct separation is not an indirect band gap. An isolated family can contain both occupied and empty states, and it can have a negative indirect gap with another energy extremum while remaining directly separated from excluded states at each k\mathbf k. Track its projector rather than energy-sorted labels through internal crossings.

For entangled source bands, Nout(k)>JN_{\rm out}(\mathbf k)>J at least somewhere. Changing rectangular VkV_{\mathbf k} changes PkP_{\mathbf k} and therefore changes the retained model. Subsequent U(J)U(J) optimization cannot undo this physical selection. Two constructions can reproduce nearly identical energies while representing different orbital content, velocities, spin matrices, or interaction projections.

This distinction determines the error ledger. Gauge-start sensitivity inside one fixed projector is an optimizer question. Sensitivity to rank, outer window, frozen window, or chemically distinct trials is representation uncertainty.

Trial functions initialize a candidate frame; they do not prove that the requested subspace exists or is unique. For trial columns GkG_{\mathbf k}, define

Ak=Ψk†Gk,Gk=Ak†Ak.A_{\mathbf k} = \Psi_{\mathbf k}^{\dagger}G_{\mathbf k}, \qquad \mathcal G_{\mathbf k} = A_{\mathbf k}^{\dagger}A_{\mathbf k}.

When AkA_{\mathbf k} has full column rank, the polar or Löwdin projection is

Vk(0)=AkGk−1/2.V_{\mathbf k}^{(0)} = A_{\mathbf k}\mathcal G_{\mathbf k}^{-1/2}.

For generalized coefficients, replace the first equation by Ak=Ck†SkGkA_{\mathbf k}=C_{\mathbf k}^{\dagger}S_{\mathbf k}G_{\mathbf k}. Report the minimum singular value of AkA_{\mathbf k} and the condition number of Gk\mathcal G_{\mathbf k} over the complete construction mesh. A chemically plausible trial set can lose rank at an uninspected momentum.

In the entangled branch, the outer window determines which source states are available and the frozen window identifies source states that must be contained in the selected subspace. Pointwise consistency requires

Nout(k)≥J,Nfr(k)≤J.N_{\rm out}(\mathbf k)\ge J, \qquad N_{\rm fr}(\mathbf k)\le J.

If FkF_{\mathbf k} contains an orthonormal frame for the frozen states in the outer coefficient space, verify

rfr(k)=∥(I−VkVk†)Fk∥2.r_{\rm fr}(\mathbf k) = \left\| \left(I-V_{\mathbf k}V_{\mathbf k}^{\dagger}\right) F_{\mathbf k} \right\|_2.

An impossible count or nonzero containment residual is a failed setup, not an optimizer failure. Widening an outer window can repair missing source states; narrowing a frozen window can repair an overconstrained target. Changing JJ changes the model and must start a new audit.

Choose trials and windows by the declared downstream object rather than by a desired chemical story. Test at least one chemically distinct admissible trial family and more than one plausible window. For spinor problems, state whether each trial is a genuine-spin orbital, a projected pseudospin, or one member of a Kramers-paired construction. Record site centers, local axes, radial functions, embedding phases, and normalization.

Energy windows alone do not encode band connectivity, symmetry content, orbital character, or operator fidelity. A narrow window may fit a few target energies while producing unstable projectors. A wide window may improve rank but add ligand or high-energy content that changes the intended low-energy model. These are scientifically meaningful alternatives, not nuisance parameters to hide.

Optimize Gauge, Spread, Symmetry, and Sewing

Section titled “Optimize Gauge, Spread, Symmetry, and Sewing”

With selected cell-periodic states, define neighbor overlaps using the physical inner product and the declared reciprocal-boundary sewing:

Mmn(k,b)=⟨umk(Φ)|un,k+b(Φ)⟩.M_{mn}^{(\mathbf k,\mathbf b)} = \left\langle u^{(\Phi)}_{m\mathbf k} \middle| u^{(\Phi)}_{n,\mathbf k+\mathbf b} \right\rangle.

The point k+b\mathbf k+\mathbf b may lie in a neighboring reciprocal cell. Its frame must first be mapped through the archived sewing convention. Treating stored components as naively periodic can create false discontinuities, incorrect centers, and broken symmetries.

For declared finite-difference weights wbw_{\mathbf b}, a standard gauge-invariant discrete spillage contribution is

ΩI=1Nk∑k,bwb[J−∑m,n=1J∣Mmn(k,b)∣2].\begin{aligned} \Omega_I &= \frac{1}{N_k} \sum_{\mathbf k,\mathbf b}w_{\mathbf b} \left[ J- \sum_{m,n=1}^{J} \left|M_{mn}^{(\mathbf k,\mathbf b)}\right|^2 \right]. \end{aligned}

For a fixed selected projector, the total quadratic spread has the form

Ω=ΩI+Ω~,\Omega = \Omega_I+\widetilde\Omega,

where Ω~\widetilde\Omega depends on the internal gauge. The entangled branch first minimizes a constrained projector functional and only afterward minimizes the gauge-dependent spread. Comparing total spreads across different projectors mixes these two tasks. A smaller spread is not automatically a more physical subspace.

Run multiple seeds and archive unsuccessful starts. Inspect the spread history, gradient or residual, individual spreads, centers, neighbor-overlap singular values, and stability under mesh refinement. A repeatable local minimum is evidence about the optimizer, not about uniqueness under windows and trials.

For a unitary symmetry represented on the source outer frame by Bg(k)B_g(\mathbf k), a symmetry-compatible frame obeys

Bg(k)Qk=QgkDg(k).B_g(\mathbf k)Q_{\mathbf k} = Q_{g\mathbf k}D_g(\mathbf k).

Antiunitary symmetries conjugate QkQ_{\mathbf k}. Nonsymmorphic translations, orbital embeddings, spin rotations, Kramers structure, and reciprocal sewing phases must be included in the residual. A symmetry-constrained solution can legitimately be less localized than an unconstrained one. The correct question is whether it meets the declared localization and symmetry tolerances, not whether it wins a spread contest against a different representation.

Transform Hamiltonians and Operators to Real Space

Section titled “Transform Hamiltonians and Operators to Real Space”

On a uniform Born–von Karman mesh, the selected frame defines

∣wnR⟩=1Nk∑ke−ik⋅R∑a=1Nout(k)∣ψak⟩Qan(k).\begin{aligned} |w_{n\mathbf R}\rangle &= \frac{1}{\sqrt{N_k}} \sum_{\mathbf k}e^{-i\mathbf k\cdot\mathbf R} \sum_{a=1}^{N_{\rm out}(\mathbf k)} |\psi_{a\mathbf k}\rangle Q_{an}(\mathbf k). \end{aligned}

This page uses that transform operationally; Wannier Functions owns its exact derivation, regularity conditions, and localization theorems.

In a source eigenbasis with diagonal energy matrix ε(k)\varepsilon(\mathbf k), set

HW(k)=Qk†ε(k)Qk,H^W(\mathbf k) = Q_{\mathbf k}^{\dagger} \varepsilon(\mathbf k) Q_{\mathbf k}, H(R)=1Nk∑ke−ik⋅RHW(k).H(\mathbf R) = \frac{1}{N_k} \sum_{\mathbf k} e^{-i\mathbf k\cdot\mathbf R} H^W(\mathbf k).

For a retained real-space set S\mathcal S, reconstruct

Hinterp(k′)=∑R∈Se+ik′⋅RH(R),H_{\rm interp}(\mathbf k') = \sum_{\mathbf R\in\mathcal S} e^{+i\mathbf k'\cdot\mathbf R}H(\mathbf R),

and verify

H(−R)=H(R)†.H(-\mathbf R)=H(\mathbf R)^{\dagger}.

For a complete isolated subspace, the untruncated transform is an exact basis change in the finite discrete problem. Continuous-momentum interpolation still carries mesh error. In an entangled branch, HWH^W is a projected Hamiltonian and need not reproduce every source eigenvalue outside the frozen or declared validation domain.

For a translation-invariant one-body operator with a valid single-momentum source matrix, use

OW(k)=Qk†Osrc(k)Qk,O^W(\mathbf k) = Q_{\mathbf k}^{\dagger} O^{\rm src}(\mathbf k) Q_{\mathbf k},

then Fourier transform with the same convention. Position, dipole, connection, two-momentum, nonlocal, current, and interaction operators can require additional derivative, embedding, augmentation, or multi-momentum terms. Transforming the Hamiltonian does not manufacture those terms.

When velocity is requested, preserve and project the physical source operator:

vαP(k)=Qk†vαsrc(k)Qk.v^P_{\alpha}(\mathbf k) = Q_{\mathbf k}^{\dagger} v^{\rm src}_{\alpha}(\mathbf k) Q_{\mathbf k}.

If the selected subspace is invariant under the represented Hamiltonian and vαsrc=ℏ−1∂kαHsrcv^{\rm src}_{\alpha}=\hbar^{-1}\partial_{k_\alpha}H^{\rm src}, define Aα=i⟨Φ∣∂kαΦ⟩\mathcal A_\alpha=i\langle\Phi|\partial_{k_\alpha}\Phi\rangle and obtain

vαP(k)=1ℏ[∂kαHW(k)−i[Aα(k),HW(k)]].\begin{aligned} v^P_{\alpha}(\mathbf k) &= \frac{1}{\hbar} \left[ \partial_{k_\alpha}H^W(\mathbf k) -i[\mathcal A_\alpha(\mathbf k),H^W(\mathbf k)] \right]. \end{aligned}

For a generic non-invariant disentangled projector, excluded-space leakage adds terms to this identity. Use and validate the source velocity matrix instead. A derivative-only effective Hamiltonian may reproduce diagonal group velocities yet miss interband, optical, spin-current, or nonlocal contributions.

Validate Full-Zone Interpolation, Truncation, and Downstream Operators

Section titled “Validate Full-Zone Interpolation, Truncation, and Downstream Operators”

Validation data must be independent of the construction mesh and must cover the full Brillouin zone. A shifted mesh is often useful because it avoids reusing the fitted points. Retain local refinement around small gaps, avoided crossings, and Fermi-surface features when those objects determine the claim.

Compare spectra as unordered sets near degeneracies. Report at least maximum and weighted RMS energy residuals over the declared validation window. A small error on a high-symmetry path does not validate a full-zone gap, pocket, projector, Fermi surface, or Brillouin-zone integral.

For a real-space cutoff RcR_c, define the auditable Hamiltonian tail

BH(Rc)=∑R∉S(Rc)∥H(R)∥2.B_H(R_c) = \sum_{\mathbf R\notin\mathcal S(R_c)} \|H(\mathbf R)\|_2.

The triangle inequality gives

∥δH(k)∥2≤BH(Rc).\|\delta H(\mathbf k)\|_2 \le B_H(R_c).

This bounds explicit truncation of the represented discrete Fourier series. It does not bound source-model error, momentum-space aliasing, subspace selection error, or omitted operator terms. Mesh density and RcR_c are coupled controls and must be converged together.

For equal-rank candidate projectors in the same physical source fiber, use

dPab(k)=∥Pa(k)−Pb(k)∥2=sin⁡θmax⁡(k).d_P^{ab}(\mathbf k) = \|P_a(\mathbf k)-P_b(\mathbf k)\|_2 = \sin\theta_{\max}(\mathbf k).

The largest principal angle detects a change of subspace that no internal U(J)U(J) gauge can remove. Report both the maximum and a distribution summary, such as the 95th percentile, so that an isolated severe failure is not hidden by an average.

Before comparing operator matrices, align the equal-rank physical frames with

Rab(k)=polar⁡[Φa†(k)Φb(k)].R_{ab}(\mathbf k) = \operatorname{polar} \left[ \Phi_a^{\dagger}(\mathbf k) \Phi_b(\mathbf k) \right].

Using the physical frames rather than their outer-space coefficients allows the two candidates to have different outer-window dimensions. The overlap must have full rank. If its smallest singular value falls below the declared correspondence tolerance, report a failed subspace match instead of forcing an operator percentage.

With normalized holdout weights ∑kwk=1\sum_{\mathbf k}w_{\mathbf k}=1, define

rOab=[∑kwk∥Oa−RabObRab†∥F2]1/2max⁡{[∑kwk∥Oa∥F2]1/2,J Ofloor}.\begin{aligned} r_O^{ab} &= \frac{ \left[ \sum_{\mathbf k}w_{\mathbf k} \left\| O_a-R_{ab}O_bR_{ab}^{\dagger} \right\|_F^2 \right]^{1/2} }{ \max\left\{ \left[ \sum_{\mathbf k}w_{\mathbf k}\|O_a\|_F^2 \right]^{1/2}, \sqrt{J}\,O_{\rm floor} \right\} }. \end{aligned}

Declare OfloorO_{\rm floor} before inspecting the comparison and archive the absolute RMS residual as well. This global normalization avoids meaningless pointwise percentages at symmetry-enforced zeros.

For a Cartesian vector operator with Nc=3N_c=3, replace each squared Frobenius norm in the numerator and denominator by the sum over α=x,y,z\alpha=x,y,z, and replace the scalar floor J Ofloor\sqrt{J}\,O_{\rm floor} by JNc Ofloor\sqrt{J N_c}\,O_{\rm floor}. All velocity and spin percentages below use this stacked Cartesian global RMS, not a componentwise maximum or a pointwise ratio.

The acceptance record should vary construction mesh, holdout mesh, real-space range, seeds, windows, trials, and rank as appropriate. Keep source-model error separate. A converged Kohn–Sham localized representation remains an auxiliary one-particle representation; numerical precision does not turn it into a quasiparticle Hamiltonian.

Diagnose Numerical Failure, Symmetry Conflict, and Topological Obstruction

Section titled “Diagnose Numerical Failure, Symmetry Conflict, and Topological Obstruction”

Diagnose failures in a fixed order.

  1. Source failure. Recheck metric positivity, electron and spinor conventions, stored band count, full-zone isolation or outer-space coverage, reciprocal sewing, source symmetries, and requested operator matrices.
  2. Setup failure. Check projection rank, Nout≥JN_{\rm out}\ge J, Nfr≤JN_{\rm fr}\le J, frozen containment, trial embedding, local axes, and prospective tolerances.
  3. Optimizer failure. Compare seeds, residual histories, gradients, neighbor-overlap conditioning, and mesh refinement while holding the projector-defining choices fixed.
  4. Representation instability. Vary admissible windows, chemically distinct trials, and rank. Compare projectors and requested operators, not only spreads and energies.
  5. Symmetry conflict. Verify source and target representations, unitary or antiunitary actions, nonsymmorphic phases, Kramers structure, orbital embeddings, and reciprocal sewing before interpreting a residual.
  6. Obstruction escalation. Only after the numerical audit passes should a failed symmetry-respecting or localized construction be routed to the existence and topology owners.

Failure to converge is not itself a topological invariant. A nonzero Chern number of a declared isolated projector gives a genuine obstruction to an exponentially localized orthonormal translated basis for precisely that projector, but Chern Numbers in Band Theory owns the invariant calculation. Time-reversal, fragile, obstructed-atomic-limit, and crystalline-symmetry obstructions impose different conditions. Symmetry of Bloch States owns irreducible representations, compatibility, and band-representation diagnostics.

Wannier centers, individual spreads, onsite energies, and hoppings are gauge, embedding, origin, or branch dependent. Berry-Phase Polarization and Charge Pumping adds filling, electron charge, ionic positions, origin, branch, and adiabatic path before turning centers into a polarization or pump claim.

Likewise, energy interpolation alone does not license interaction or response physics. Hubbard Physics in Materials owns active-space adequacy, screened interactions, double counting, and correlated-model validation. Transport owners retain lifetimes, collision integrals, vertex corrections, contacts, and dc orders of limits.

This synthetic record inherits the source problem from Band Structure Workflows. It audits a representation; it is not a prediction for a real compound.

  1. Physical problem. Use the fixed tetragonal P4/mmmP4/mmm (No. 123) Bi2Te2\mathrm{Bi_2Te_2} teaching crystal with a=4.10 A˚a=4.10\ \text{Å}, c=6.20 A˚c=6.20\ \text{Å}, Bi at (0,0,0.22)(0,0,0.22) and (0,0,0.78)(0,0,0.78), and Te at (1/2,1/2,0.28)(1/2,1/2,0.28) and (1/2,1/2,0.72)(1/2,1/2,0.72). The target is the infinite periodic bulk primitive cell with its archived origin and embedding.
  2. State and limit. Use the neutral fixed-NN, nonmagnetic, inversion- and time-reversal-symmetric PBE+SOC source at zero physical temperature, field, strain, and drive. The fixed structure contains 62 explicitly enumerated spinor electrons and uses gs=1g_s=1.
  3. Claim and accuracy. Construct a symmetry-compatible localized representation of the exported rank-88 near-gap spinor manifold. Require maximum holdout energy error below 2 meV2\ \text{meV}, global projected- velocity RMS below 1%1\%, and BH<1 meVB_H<1\ \text{meV}. No material, fundamental-gap, optical, topological, polarization, or transport claim is requested.
  4. Representation and provenance. Consume immutable artifact BITE-P4MMM-PBE-SOC-062E-096B-v1. It contains the accepted PBE+SOC density, 96 stored spinor bands, eigenvalues and overlaps, unit-cell and sewing data, symmetry matrices, and physical velocity matrices including the nonlocal- pseudopotential contribution. Archive the teaching identifiers, code and environment versions, and pseudopotential hashes.
  5. Method and controlled domain. The eight-band source projector has minimum direct separation Δsep=0.36 eV\Delta_{\rm sep}=0.36\ \text{eV} from excluded bands at every momentum. Set Nout=J=8N_{\rm out}=J=8. Outer and frozen windows, frozen containment, and disentanglement are not applicable. Start from four site-centered pzp_z-like spatial trials and their Kramers partners, use a full-rank polar projection, and perform symmetry-constrained spread minimization with six seeds.
  6. Finite numerical problem. Use unshifted Γ\Gamma-centered construction meshes 12312^3, 16316^3, and 20320^3. The independent holdout contains every point kn=∑i(ni+1/2)bi/23\mathbf k_{\mathbf n}=\sum_i(n_i+1/2)\mathbf b_i/23 with ni=0,…,22n_i=0,\ldots,22, each with weight 23−323^{-3}. Use real-space cutoffs Rc=12,15,18 A˚R_c=12,15,18\ \text{Å}. Store every source state, neighbor link, seed, and discarded real-space matrix.
  7. Estimator and forward model. Record total and individual spreads, centers, minimum trial singular value, maximum and RMS holdout spectral errors, gauge-aligned global velocity RMS, BHB_H, Hermiticity, symmetry and Kramers sewing, and artifact hashes. Declare vfloor=10−3Espana/ℏv_{\rm floor}=10^{-3}E_{\rm span}a/\hbar with Espan=2.0 eVE_{\rm span}=2.0\ \text{eV} and retain the absolute velocity RMS as well as its percentage. Compare spectra as sets and operators as subspaces near degeneracies.
  8. Convergence and uncertainty. The three meshes give Ω=24.8,24.1,24.0 A˚2\Omega=24.8,24.1,24.0\ \text{Å}^2; maximum holdout energy residuals 8.0,3.1,1.4 meV8.0,3.1,1.4\ \text{meV}; RMS energy residuals 2.4,0.9,0.38 meV2.4,0.9,0.38\ \text{meV}; global velocity RMS residuals 2.6%,1.1%,0.48%2.6\%,1.1\%,0.48\%; and minimum trial singular values 0.44,0.46,0.460.44,0.46,0.46. At the accepted mesh, BH(12,15,18 A˚)=5.8,2.1,0.76 meVB_H(12,15,18\ \text{Å})=5.8,2.1,0.76\ \text{meV}. Final center changes are below 0.002 A˚0.002\ \text{Å} and individual-spread changes below 0.02 A˚20.02\ \text{Å}^2. Source-functional uncertainty remains separate from these representation errors.
  9. Verification, validation, and provenance. The accepted run has Hermiticity residual below 10−1210^{-12}, passed orthonormality and unitarity checks, symmetry-sewing residual below 10−810^{-8}, correct time-reversal and Kramers sewing, and independent holdout evidence. Source matrices, seeds, failed starts, scripts, versions, hashes, and the absolute operator residual are archived.
  10. Licensed claim and stopping rule. License only a localized representation of the declared auxiliary Kohn–Sham rank-88 subspace at the stated energy, velocity, and tail tolerances. Loss of direct isolation, projection rank, sewing, holdout, or operator accuracy stops the export. Invariant, polarization, transport, and correlated-model conclusions go to their canonical owners.

All three acceptance tolerances pass only at the 20320^3 construction mesh and 18 A˚18\ \text{Å} cutoff. The result is useful precisely because its license is narrow: it validates one auxiliary representation and one requested operator, not every physical interpretation that could be attached to the files.

Worked Audit: An Entangled Metallic Manifold

Section titled “Worked Audit: An Entangled Metallic Manifold”

This second synthetic record separates numerical convergence from stability of the selected representation.

  1. Physical problem. Use a hypothetical bulk tetragonal P4/mmmP4/mmm (No. 123) MX2MX_2 teaching crystal with a=3.80 A˚a=3.80\ \text{Å}, c=6.00 A˚c=6.00\ \text{Å}, MM at (0,0,0)(0,0,0), and XX at (1/2,0,1/2)(1/2,0,1/2) and (0,1/2,1/2)(0,1/2,1/2).
  2. State and limit. Use a neutral fixed-NN, nonmagnetic, inversion- and time-reversal-symmetric infinite periodic PBE+SOC source at zero physical temperature, field, strain, and drive. The archived chemical potential is the zero of energy. The source-data owner has completed and stored the metallic occupation and smearing-to-zero audit.
  3. Claim and accuracy. Test a rank-66 t2gt_{2g}-like spinor representation crossing the chemical potential. Require maximum frozen-window energy error at most 5 meV5\ \text{meV}, maximum projector distance below 0.050.05 under admissible window and trial changes, and global velocity and spin RMS changes below 2%2\%.
  4. Representation and provenance. Consume immutable source artifact MX2-P4MMM-SOC-028E-064B-v1, with 28 explicit spinor electrons per cell in 64 stored spinor bands. It records fully relativistic norm-conserving PBE+SOC teaching pseudopotentials with 12-valence MM and 8-valence XX, a 90 Ry90\ \text{Ry} cutoff, the accepted 18318^3 source density, overlaps and sewing matrices, physical velocity including nonlocal terms, and S=ℏσ/2\mathbf S=\hbar\boldsymbol\sigma/2 matrices. Archive all source hashes. Wannierization never updates the density.
  5. Method and controlled domain. Select and localize a rank-66 projector with frozen-state containment. Begin with three t2gt_{2g} spatial trials and their Kramers partners, compare them with d+pd+p-tailed trials, and use twelve starts. The output is only a candidate one-body representation, not a correlated, transport, or topological model.
  6. Finite numerical problem. Use unshifted Γ\Gamma construction meshes 10310^3, 14314^3, and 18318^3. The independent holdout contains every point kn=∑i(ni+1/2)bi/23\mathbf k_{\mathbf n}=\sum_i(n_i+1/2)\mathbf b_i/23 with ni=0,…,22n_i=0,\ldots,22, each with weight 23−323^{-3}. Use Rc=12,15,18 A˚R_c=12,15,18\ \text{Å}. The baseline outer and frozen windows are [−3,2][-3,2] and [−0.6,0.6] eV[-0.6,0.6]\ \text{eV} relative to the chemical potential. Test outer windows [−2.5,1.5][-2.5,1.5] and [−3.5,2.5] eV[-3.5,2.5]\ \text{eV}, frozen windows [−0.4,0.4][-0.4,0.4] and [−0.8,0.8] eV[-0.8,0.8]\ \text{eV}, and both trial families. Across the mesh, 8≤Nout(k)≤148\le N_{\rm out}(\mathbf k)\le14 and Nfr(k)≤6N_{\rm fr}(\mathbf k)\le6.
  7. Estimator and forward model. Record spread, maximum and RMS holdout energy residuals, BHB_H, frozen containment, Hermiticity, symmetry sewing, maximum and 95th-percentile dPd_P, and gauge-aligned global operator RMS. Use vfloor=10−3Espana/ℏv_{\rm floor}=10^{-3}E_{\rm span}a/\hbar with Espan=5 eVE_{\rm span}=5\ \text{eV} and Sfloor=10−3ℏ/2S_{\rm floor}=10^{-3}\hbar/2. Archive absolute RMS values and never form pointwise percentages at symmetry-enforced zeros.
  8. Convergence and uncertainty. The baseline meshes give Ω=19.6,19.0,18.9 A˚2\Omega=19.6,19.0,18.9\ \text{Å}^2, maximum holdout energy residuals 13,6,3 meV13,6,3\ \text{meV}, and RMS energy residuals 4.1,1.8,0.9 meV4.1,1.8,0.9\ \text{meV}. At the accepted mesh, BH(12,15,18 A˚)=6.4,2.3,0.82 meVB_H(12,15,18\ \text{Å})=6.4,2.3,0.82\ \text{meV}, frozen containment is below 10−1010^{-10}, symmetry sewing below 10−810^{-8}, and Hermiticity below 10−1210^{-12}. Nevertheless, baseline and wide admissible windows differ by maximum dP=0.42d_P=0.42 with 95th percentile 0.060.06. Alternate trials give Ω=18.8\Omega=18.8 versus 18.9 A˚218.9\ \text{Å}^2 but maximum dP=0.38d_P=0.38. Gauge-aligned velocity and spin matrices change by global RMS 8%8\% and 12%12\%, and gauge-qualified onsite-matrix eigenvalues shift by 0.11 eV0.11\ \text{eV}. The onsite shift is an ancillary basis-dependent diagnostic, not an invariant.
  9. Verification, validation, and provenance. Optimizer, frozen containment, mesh, tail, Hermiticity, and symmetry checks pass, but representation validation fails. Archive all alternatives rather than only the smallest-spread result. A diagnostic rank-88 rerun with augmented ligand character gives maximum dP=0.04d_P=0.04, maximum velocity or spin global RMS change 1.5%1.5\%, maximum energy residual 2 meV2\ \text{meV}, and tail below 1 meV1\ \text{meV}. Because rank changed, this is a new candidate and still requires its own complete ten-field audit.
  10. Licensed claim and stopping rule. Do not export the rank-66 representation for Hubbard, topology, or transport work even though its energy interpolation converged. Omitted ligand weight is the credible failure mode. The rank-88 representation may be licensed only after its fresh audit confirms every declared tolerance.

This audit is a deliberate non-result. Its energy, optimizer, tail, and symmetry tests pass, yet the selected rank-66 physical subspace and its operators are not stable. Reporting the attractive energy overlay would hide the failure that matters downstream.

Exit Checkpoint, Artifacts, and Canonical Handoffs

Section titled “Exit Checkpoint, Artifacts, and Canonical Handoffs”

Before exporting a localized representation, verify:

  • the source physical problem and artifact identifiers remained fixed;
  • isolated or entangled status was established over the full zone;
  • rank, windows, frozen constraints, trials, spinor conventions, and symmetry requirements were declared before optimization;
  • projection rank and conditioning passed on the complete construction mesh;
  • multiple seeds and chemically distinct admissible trials were retained;
  • energies, projectors, real-space tails, sewing, and every requested operator passed prospective tolerances on independent full-zone data;
  • mesh and real-space range were converged together;
  • basis-dependent centers, spreads, onsite terms, and hoppings were labelled as such;
  • failed runs and credible alternative representations remain in the archive; and
  • the final sentence licenses only the tested representation and downstream use.

The reproducibility package should include source and output identifiers, input structures, cell and reciprocal conventions, raw states and overlaps, operator matrices, windows, trials, seeds, symmetry settings, convergence histories, rejected alternatives, real-space matrices, holdout data, analysis scripts, code and environment versions, and cryptographic hashes. Preserve raw matrices rather than only plots or formatted summaries.

Canonical handoffs keep the computational result from becoming an overclaim:

Reusable algorithms, package engineering, performance, production notebooks, and benchmarks remain with Computational QM. This page records enough mathematics to audit a workflow; it does not duplicate those future implementation owners.

Classify each step as an exact basis change, a definition of a reduced subspace, a representation choice, a numerical approximation, or a physical modeling approximation: a U(J)U(J) rotation inside a fixed projector; the complete finite Bloch–Wannier transform; projection onto an isolated band family; entangled VkV_{\mathbf k} selection; finite momentum sampling; real-space truncation; and projection of Coulomb interactions.

Solution

A U(J)U(J) rotation and the complete finite transform are exact basis changes inside a fixed selected subspace. Projection onto an isolated band family defines the retained spectral projector exactly, but discarding its complement is a physical reduction whose adequacy depends on the requested scale and observable.

Entangled VkV_{\mathbf k} selection changes the projector and is a representation or model choice. Finite sampling and real-space truncation are numerical approximations to the continuous and untruncated representation. Projecting Coulomb interactions is exact only if every transformed matrix element and all eliminated-sector effects are retained. Practical screened few-parameter interactions add physical modeling, screening, and double-counting choices owned by Hubbard Physics in Materials.

At one momentum, a three-state outer space and two trial functions give

A=(1000.2000).A= \begin{pmatrix} 1&0\\ 0&0.20\\ 0&0 \end{pmatrix}.

Compute the singular values, G=A†A\mathcal G=A^{\dagger}A, its 2-norm condition number, and V(0)=AG−1/2V^{(0)}=A\mathcal G^{-1/2}. Is the polar projection defined? A prospective workflow requires σmin⁡≥0.10\sigma_{\min}\ge0.10 everywhere. What happens if mesh refinement lowers the second singular value to 0.040.04?

Solution

The singular values are 11 and 0.200.20, and

G=(1000.04),κ2(G)=25.\mathcal G = \begin{pmatrix} 1&0\\ 0&0.04 \end{pmatrix}, \qquad \kappa_2(\mathcal G)=25.

Both singular values are nonzero, so the polar projection is mathematically defined:

V(0)=(100100).V^{(0)} = \begin{pmatrix} 1&0\\ 0&1\\ 0&0 \end{pmatrix}.

The initial point also passes the declared 0.100.10 threshold, although the second direction is much less robust. A refined-mesh singular value 0.040.04 fails the prospective requirement. The trial set does not span the target reliably over the full mesh; keeping the coarse result or adding optimizer iterations is not a repair. Change the trials or the declared target subspace and repeat the audit.

A proposed entangled calculation uses J=6J=6. At one momentum the outer window contains five states. At another, the frozen window contains seven. Diagnose both failures and give valid repairs that do not disguise a rank change as ordinary convergence.

Solution

The first point violates Nout≥JN_{\rm out}\ge J because 5<65<6; no rank-six semiunitary selection can be built from five source states. Widen the outer window or supply more source bands while holding the source physical problem fixed.

The second point violates Nfr≤JN_{\rm fr}\le J because seven independent frozen states cannot fit inside a rank-six projector. Narrow the frozen window or revise which states are mandatory. Increasing JJ is also possible, but it defines a new representation and requires a fresh ten-field audit. Neither failure is cured by changing the spread-minimization tolerance or running more seeds.

Exercise 4: Reconstruct and truncate a Hamiltonian

Section titled “Exercise 4: Reconstruct and truncate a Hamiltonian”

On a one-dimensional four-point mesh k=0,π/2,π,3π/2k=0,\pi/2,\pi,3\pi/2, a one-orbital Wannier-frame Hamiltonian has values HW(k)=2,0,−2,0 eVH^W(k)=2,0,-2,0\ \text{eV}. Using

H(R)=14∑ke−ikRHW(k),H(R)=\frac14\sum_k e^{-ikR}H^W(k),

find the independent coefficients H(0)H(0), H(1)H(1), H(−1)H(-1), and H(2)H(2). Verify Hermiticity. If a faulty truncation keeps only R=0,1R=0,1, compute its BHB_H and explain the defect.

Solution

Direct evaluation gives

H(0)=0,H(1)=H(−1)=1 eV,H(2)=0.H(0)=0, \qquad H(1)=H(-1)=1\ \text{eV}, \qquad H(2)=0.

The coefficients satisfy H(−R)=H(R)∗H(-R)=H(R)^*, as required for this scalar problem. Keeping R=1R=1 while discarding its Hermitian partner R=−1R=-1 gives

BH=∣H(−1)∣=1 eV.B_H=|H(-1)|=1\ \text{eV}.

The reconstructed function is then complex and the truncation violates Hermiticity. A symmetry-respecting cutoff keeps both R=±1R=\pm1; for the supplied finite Fourier series its remaining tail is zero.

Three construction meshes give maximum energy errors 0.6,0.4,0.3 meV0.6,0.4,0.3\ \text{meV} on a conventional path. On an independent shifted full-zone mesh, the corresponding maximum errors are 12.0,5.2,2.4 meV12.0,5.2,2.4\ \text{meV} and RMS errors are 3.8,1.4,0.7 meV3.8,1.4,0.7\ \text{meV}. The prospective maximum-error tolerance is 2 meV2\ \text{meV}. Which mesh is accepted, and what should happen next?

Solution

None is accepted. Every path result looks excellent, but the full-zone maximum remains above 2 meV2\ \text{meV} even on the finest mesh. The decreasing RMS suggests improvement, yet an off-path avoided crossing or extremum still controls the maximum claim.

Refine the construction mesh and the offending full-zone region while keeping the source Hamiltonian fixed. Recheck the selected projector near the large residual and repeat the independent holdout. The path cannot be promoted to the validation estimator because it never samples the failed region.

At one momentum in an invariant two-band subspace, suppose

HW=(0001)eV,H^W= \begin{pmatrix} 0&0\\ 0&1 \end{pmatrix}\text{eV}, ∂kHW=(20.30.3−1)eV A˚,A=(00.2i−0.2i0)A˚.\partial_k H^W= \begin{pmatrix} 2&0.3\\ 0.3&-1 \end{pmatrix}\text{eV Å}, \qquad \mathcal A= \begin{pmatrix} 0&0.2i\\ -0.2i&0 \end{pmatrix}\text{Å}.

Compute ℏv=∂kHW−i[A,HW]\hbar v=\partial_kH^W-i[\mathcal A,H^W]. Which entries are group velocities, and why is this calculation insufficient for a generic disentangled projector?

Solution

The commutator contribution is

−i[A,HW]=(00.20.20)eV A˚,-i[\mathcal A,H^W] = \begin{pmatrix} 0&0.2\\ 0.2&0 \end{pmatrix}\text{eV Å},

so

ℏv=(20.50.5−1)eV A˚.\hbar v = \begin{pmatrix} 2&0.5\\ 0.5&-1 \end{pmatrix}\text{eV Å}.

In the displayed eigenbasis of HWH^W, the diagonal entries divided by ℏ\hbar are the two group velocities. The off-diagonal 0.5 eV A˚/ℏ0.5\ \text{eV Å}/\hbar is an interband matrix element and would be lost by differentiating eigenvalues alone.

The identity used the assumption that the selected subspace is invariant under the represented Hamiltonian. A generic disentangled projector leaks into its excluded complement, adding terms not recoverable from HWH^W and its internal connection. In that case project and validate the archived physical source velocity, including nonlocal or augmentation contributions.

Two equal-rank candidate frames on the same holdout point have physical-frame overlap singular values 0.9990.999, 0.9800.980, and 0.7000.700. Compute the largest projector distance. Can an internal gauge rotation repair the difference if the prospective tolerance is dP<0.05d_P<0.05?

Solution

The smallest overlap singular value is the cosine of the largest principal angle. Therefore

dP=sin⁡θmax⁡=1−0.7002≃0.714.d_P = \sin\theta_{\max} = \sqrt{1-0.700^2} \simeq 0.714.

This greatly exceeds 0.050.05. An internal unitary rotation changes frames inside their projectors but leaves every principal angle between the two subspaces unchanged. It cannot repair the mismatch. Even identical energy interpolation would not make these candidates the same physical subspace; the window or trial dependence must be resolved or the export stopped.

Exercise 8: Complete a downstream export ledger

Section titled “Exercise 8: Complete a downstream export ledger”

A group proposes to use the rank-66 entangled audit above to publish a DOS, Fermi surface, dc Hall coefficient, and three-orbital Hubbard model. Complete the ten-field record, state what this page licenses, and route every physical claim.

Solution
  1. Physical problem. The declared object is the infinite periodic tetragonal MX2MX_2 teaching crystal with its fixed cell, origin, orbital embedding, and bulk boundary conditions. No surface, disorder realization, contact geometry, or experimental specimen has been supplied.
  2. State and limit. The source is neutral, fixed-NN, nonmagnetic, inversion- and time-reversal-symmetric PBE+SOC at zero physical temperature, field, strain, and drive with archived μ=0\mu=0. A dc transport order of limits and scattering state are unresolved.
  3. Claim and accuracy. The proposed DOS, Fermi surface, Hall coefficient, and Hubbard model are four different claims. Only the representation tolerances were declared. Probe resolution, DOS normalization, Fermi-sheet accuracy, Hall units and limits, and correlated-observable tolerances are unresolved.
  4. Representation and provenance. The source artifact is MX2-P4MMM-SOC-028E-064B-v1; the candidate is the rank-66 disentangled subspace with its recorded windows, trials, metric, spinor convention, velocity and spin matrices, source hashes, and omitted ligand weight.
  5. Method and controlled domain. The one-body PBE+SOC source and disentanglement workflow can propose an auxiliary localized representation. They do not supply lifetimes, collision or vertex physics, screened interactions, double counting, or a many-body solver.
  6. Finite numerical problem. The record includes 10310^3, 14314^3, and 18318^3 construction meshes and the full shifted holdout kn=∑i(ni+1/2)bi/23\mathbf k_{\mathbf n}=\sum_i(n_i+1/2)\mathbf b_i/23 for ni=0,…,22n_i=0,\ldots,22, with weights 23−323^{-3}. It also includes three real-space cutoffs, the window variants, two trial families, twelve seeds, and 8≤Nout≤148\le N_{\rm out}\le14 with Nfr≤6N_{\rm fr}\le6.
  7. Estimator and forward model. Energy, tail, projector, velocity, spin, containment, and symmetry metrics were evaluated. DOS quadrature, Fermi-surface reconstruction, scattering and Hall kernels, interaction projection, screening, and probe forward models remain unresolved.
  8. Convergence and uncertainty. Energy, mesh, tail, optimizer, containment, and symmetry controls pass, but admissible windows and trials change the projector by 0.420.42 or 0.380.38, velocity by 8%8\%, spin by 12%12\%, and onsite eigenvalues by 0.11 eV0.11\ \text{eV}. These representation errors dominate the requested downstream uses.
  9. Verification, validation, and provenance. Raw alternatives and hashes are archived, but the prospective projector and operator tolerances fail. The diagnostic rank-88 calculation is a changed model and is not a repair of the rank-66 record.
  10. Licensed claim and stopping rule. This workflow licenses no rank-66 downstream export. Stop. After a fresh accepted representation audit, route DOS normalization to Density of States, sheet geometry to Fermi Surface, dc and Hall limits to the Transport, Response, and Optics gateway, and active-space, interaction, screening, and double-counting claims to Hubbard Physics in Materials.
  • N. Marzari and D. Vanderbilt, “Maximally Localized Generalized Wannier Functions for Composite Energy Bands,” Physical Review B 56, 12847–12865 (1997), doi:10.1103/PhysRevB.56.12847.
  • N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, “Maximally Localized Wannier Functions: Theory and Applications,” Reviews of Modern Physics 84, 1419–1475 (2012), doi:10.1103/RevModPhys.84.1419.
  • A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, “wannier90: A Tool for Obtaining Maximally-Localised Wannier Functions,” Computer Physics Communications 178, 685–699 (2008), doi:10.1016/j.cpc.2007.11.016.
  • G. Pizzi et al., “Wannier90 as a Community Code: New Features and Applications,” Journal of Physics: Condensed Matter 32, 165902 (2020), doi:10.1088/1361-648X/ab51ff.
  • R. Sakuma, “Symmetry-Adapted Wannier Functions in the Maximal Localization Procedure,” Physical Review B 87, 235109 (2013), doi:10.1103/PhysRevB.87.235109.
  • I. Souza, N. Marzari, and D. Vanderbilt, “Maximally Localized Wannier Functions for Entangled Energy Bands,” Physical Review B 65, 035109 (2001), doi:10.1103/PhysRevB.65.035109.
  • X. Wang, J. R. Yates, I. Souza, and D. Vanderbilt, “Ab Initio Calculation of the Anomalous Hall Conductivity by Wannier Interpolation,” Physical Review B 74, 195118 (2006), doi:10.1103/PhysRevB.74.195118; erratum, Physical Review B 76, 169902 (2007), doi:10.1103/PhysRevB.76.169902.
  • J. R. Yates, X. Wang, D. Vanderbilt, and I. Souza, “Spectral and Fermi Surface Properties from Wannier Interpolation,” Physical Review B 75, 195121 (2007), doi:10.1103/PhysRevB.75.195121.