Computational Quantum Matter

Computational quantum matter connects models of electrons, spins and lattices to quantities such as energies, correlations and spectra. A useful calculation answers a specific physical question and shows how its answer changes when numerical settings or physical approximations change. This chapter introduces three materials calculations: band structure, Wannier interpolation and dynamical mean-field theory (DMFT).

Required background. Choosing a Model explains how to select degrees of freedom, energy scales and approximations. Here you will use those choices to decide what to calculate.

Helpful background. Convergence Tests explains refinement and extrapolation. How Quantum Matter Is Measured helps connect a calculated quantity with an experimental signal.

Choose a starting calculation

Start with the quantity you want to understand. For example, “Does this Hamiltonian have an insulating band filling?” is more precise than “Calculate the bands.” It specifies a question that needs the occupied and empty states throughout the Brillouin zone, not just a visually attractive path plot.

Your question	Start here	What you will learn to check
Bands of a crystal	Band Structure Workflows	Structure, basis, mesh and band interpretation.
A compact orbital model	Wannierization Workflows	Energy windows, interpolation and operators.
Local dynamical correlations	DMFT for Quantum Materials	Correlated orbitals, interactions, impurity solutions and spectra.

These are starting points for different questions, not a compulsory three-step sequence. A tight-binding model can be studied without a preceding density-functional calculation. Wannier functions can provide a compact basis and efficient interpolation (Marzari et al., 2012), but changing representation does not by itself add missing interactions. Material DMFT introduces a frequency-dependent local self-energy; its treatment of local correlations does not make spatial correlations exact (Kotliar et al., 2006).

For a spin model or interacting lattice Hamiltonian, use Computational Many-Body QM to compare exact diagonalization, Monte Carlo and tensor networks. The useful method depends on dimension, temperature, entanglement, sign structure and the observable. The name of a method alone does not establish its accuracy in a new regime.

Separate the physical quantity from the numerical output

An eigenvalue, an expectation value and a measured intensity answer different questions. A Kohn–Sham eigenvalue belongs to an auxiliary one-electron problem. An electron-addition or electron-removal spectrum describes excitations of an interacting system. Photoemission adds transition matrix elements, occupation factors and instrumental resolution. A resemblance between two curves is useful, but their relationship must come from the physics of the measurement. See Spectral Functions and Angle-Resolved Photoemission.

Likewise, a finite cluster has a discrete spectrum; a bulk material may have a continuum. A terminal conductance includes leads and contacts; a bulk conductivity describes a different geometry. Transport, Response, and Optics explains which response and limits fit each experiment.

Before a long run, write down the system, state and target quantity in ordinary physical terms. For example: “Find the zero-temperature ground-state energy per site of this spin chain, and extrapolate its length dependence.” This immediately suggests checks of the ground-state search, finite-size dependence and energy normalization. It also tells you what the result will not yet answer, such as the finite-temperature susceptibility.

Numerical convergence and physical accuracy

Three questions help interpret a computed number:

1. Was the finite problem solved accurately? Tighten the eigensolver or self-consistency tolerance; check residuals, normalization and symmetries.
2. Does the finite problem approximate the intended model? Enlarge the basis or system, refine the mesh, or reduce the time step. Check the actual observable of interest, since different quantities converge at different rates.
3. Is that model adequate for the physical question? Test the relevant interactions, magnetic states, temperature, disorder and other assumptions.

A smaller solver tolerance addresses the first question. It cannot repair an omitted orbital, a small simulation cell or an unsuitable interaction model. Agreement between independent codes is valuable for finding implementation errors. The DFT comparison by Lejaeghere et al. (2016), for example, compared equations of state computed using the same exchange-correlation approximation; such agreement does not remove the approximation’s physical limitations.

For a scalar observable, this distinction can be written as

$$Q_{\rm phys}-\widehat Q ={}(Q_{\rm phys}-Q_{\rm model}) +(Q_{\rm model}-Q_{\rm finite}) +(Q_{\rm finite}-\widehat Q).$$

Here $Q_{\rm finite}$ is the exact answer to the chosen finite approximation, and $\widehat Q$ is its computed estimate. The terms represent physical-model, finite-representation and numerical-solution discrepancies. The identity separates these sources of disagreement without supplying their sizes or distributions. The terms need not be independent random errors. Add variances in quadrature only when their covariance can be neglected for a justified reason.

Refine controls separately before combining them. For instance, converge a momentum mesh at several basis cutoffs; two simultaneous changes can accidentally cancel. A stable total energy also need not imply stable small energy differences, forces or response functions. Conditioning and Stability explains why sensitive quantities may need tighter control.

A perfectly solved matrix can omit an important orbital

Consider two orthonormal orbitals with the real symmetric Hamiltonian

$$H=\begin{pmatrix}0&v\\v&\Delta\end{pmatrix}, \qquad \Delta>0.$$

The coupling $v$ mixes the lower bare orbital with an orbital at energy $\Delta$. If you retain only the first orbital and simply project $H$, the matrix is $(0)$. Its energy is exactly zero: no increase in arithmetic precision or solver iterations can change that answer.

The two-orbital ground energy is instead

$$E_- = \frac{\Delta-\sqrt{\Delta^2+4v^2}}{2}.$$

For $|v|\ll\Delta$, the omitted orbital lowers the energy by approximately $v^2/\Delta$, so $E_-\simeq-v^2/\Delta$. With $\Delta=1\ \mathrm{eV}$ and $v=0.3\ \mathrm{eV}$, the full result is $-0.0831\ \mathrm{eV}$; the second-order estimate is $-0.0900\ \mathrm{eV}$. Both differ substantially from the projected answer zero. This example separates a basis omission from an approximation to the effect of that omission.

The residual makes the distinction precise. For a normalized approximate eigenvector $\widetilde\psi$, define

$$r=H\widetilde\psi-\widetilde E\widetilde\psi.$$

The vector $(1,0)^T$ with energy zero has residual zero in the projected one-orbital problem, but residual $(0,v)^T$ in the full problem. Always state which Hamiltonian the residual tests. For a finite Hermitian $H$, the spectral theorem gives

$$\min_j|E_j-\widetilde E|\leq\lVert r\rVert_2.$$

This bounds the distance to some eigenvalue of that matrix. It does not identify the ground state: an exact excited eigenvector also has zero residual. Eigenvector accuracy additionally depends on separation from other eigenvalues; near degeneracy, a subspace may be better determined than a particular vector. See Higham (2002) for numerical error analysis.

Interpreting statistical and self-consistent calculations

Monte Carlo samples are usually correlated. After equilibration, let $\rho(t)$ be the normalized autocorrelation of the measured observable at lag $t$ and set

$$\tau_{\rm int}=\frac12+\sum_{t=1}^{\infty}\rho(t).$$

For a stationary record of $N$ samples with finite variance and summable correlations, the large-$N$ variance of its mean is approximately

$$\operatorname{Var}(\overline O) \simeq\frac{2\tau_{\rm int}}{N}\operatorname{Var}(O).$$

Thus correlated measurements can still give the correct mean, but counting them as independent understates uncertainty when correlations are positive. The correlation time depends on the observable. Long runs and binning checks are needed to resolve slow modes (Sandvik, 2010); a sign-reweighted ratio also requires the numerator–denominator covariance. More samples do not correct a chain trapped in the wrong part of its state space.

Self-consistency has another limitation. A small value of $\lVert F(x)-x\rVert$ means that an iterate nearly solves the chosen fixed-point equation. It does not establish uniqueness, physical stability or convergence of every other approximation. Try different initial states and compare the appropriate thermodynamic potentials when they are available. Similarly, tensor-network discarded weight is a useful truncation diagnostic, not a universal bound on every observable’s error.

When comparing with data, preserve enough information to repeat the calculation: the model parameters and input structure, code version and numerical settings, sampling data when relevant, and the script that produces the plotted quantity. This lets a surprising feature be checked without reconstructing the entire calculation from a figure. Labs provides computational investigations; Data Interpretation and Pitfalls discusses competing physical explanations.

Exercises

1. Recover the effect of the second orbital

For the two-orbital example, eliminate the second component of an eigenvector. Derive an energy-dependent equation for the first component. How does the second-order energy estimate follow?

Solution

Writing the vector as $(a,b)^T$, the second row gives $b=-va/(\Delta-E)$ when $E\ne\Delta$. Substituting into the first row gives

$$E=-\frac{v^2}{\Delta-E}.$$

The low-energy root is $E_-$. For $|v|/\Delta\ll1$, its magnitude is much smaller than $\Delta$, so replacing $\Delta-E$ by $\Delta$ gives $E\simeq-v^2/\Delta$. Simply deleting the second orbital omits this correction. For $v=0.3\ \mathrm{eV}$ and $\Delta=1\ \mathrm{eV}$, the second-order answer differs from the exact answer by about $6.9\ \mathrm{meV}$.

2. Diagnose apparent convergence

A dimensionless illustrative estimate depends on a mesh spacing $h$ and a positive basis parameter $b$ as $Q(h,b)=2+h^2-b^{-2}$. A sequence uses $h=1/b$ for $b=10,20,40$ and always returns $2$. Does this establish control of both errors? What additional runs would help?

Solution

The two finite errors cancel on that sequence. The limiting value happens to be correct, but the sequence alone does not reveal either contribution. Hold $b=10$ and use $h=0.10,0.05,0.025$: the estimates are $2,1.9925,1.990625$. Hold $h=0.10$ and use $b=10,20,40$: they are $2,2.0075,2.009375$. Independent refinement exposes the two error directions. In a real calculation the error formula is usually unknown, which makes these separate checks more useful than the apparent agreement on one sequence.

3. Count effective samples

A stationary record has $N=200{,}000$ samples and $\tau_{\rm int}=40$ in the convention above. Estimate its effective independent sample count and the factor by which an independence assumption underestimates the standard error. Name one condition that would invalidate the estimate.

Solution

$$N_{\rm eff}\simeq\frac{N}{2\tau_{\rm int}}=2500.$$

The standard error is larger by $\sqrt{80}\simeq8.94$. This assumes the record is long enough to resolve the relevant correlations and the chain samples the desired stationary distribution. An unresolved slow transition, an unequilibrated chain or a ratio estimator treated without its covariance can invalidate it.

4. Separate a band separation from an insulating gap

Consider a one-dimensional two-band model $E_\pm(k)=2t\cos k\pm\Delta/2$, with $t,\Delta>0$ and dimensionless $k\in[-\pi,\pi]$. At filling one spinless electron per cell, is a positive same-momentum separation $\Delta$ sufficient for a band insulator? Evaluate the indirect gap when $\Delta=2t$.

Solution

The bottom of the upper band is $-2t+\Delta/2$ and the top of the lower band is $2t-\Delta/2$. Their difference is

$$E_{\rm g}^{\rm ind}=\Delta-4t.$$

For $\Delta=2t$ it is $-2t$: the bands overlap in energy, despite being separated at each momentum. The lowest-energy filling occupies parts of both bands, so the model is not a band insulator. A positive indirect gap requires $\Delta>4t$. This is a statement about the specified noninteracting model; connecting a Kohn–Sham calculation to an excitation gap requires further physical reasoning, as discussed in Band Structure Workflows.

References

• N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002.
• G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, “Electronic Structure Calculations with Dynamical Mean-Field Theory,” Reviews of Modern Physics 78, 865–951 (2006), doi:10.1103/RevModPhys.78.865.
• K. Lejaeghere et al., “Reproducibility in Density Functional Theory Calculations of Solids,” Science 351, aad3000 (2016), doi:10.1126/science.aad3000.
• 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. W. Sandvik, “Computational Studies of Quantum Spin Systems,” AIP Conference Proceedings 1297, 135–338 (2010), doi:10.1063/1.3518900, author’s notes.