Questaal: a package of electronic structure methods based on the linear muffin-tin orbital technique

PROGRAM SUMMARY

• Questaal calculates electronic structure and resulting physical, spectroscopic, and magnetic properties for periodic solids across diverse material classes and kinds and strengths of electronic correlation.Its stated nature is highly accurate ab initio calculation.Source: PROGRAM SUMMARY · p. 1
• Questaal is licensed under GNU General Public License version 3 and is written in Fortran, C, Python, and Shell.Source: PROGRAM SUMMARY · p. 1
• The many-electron problem is treated using density functional theory, GW many-body perturbation theory at varying self-consistency levels, notably quasiparticle self-consistent GW, and dynamical mean field theory.Source: PROGRAM SUMMARY · p. 1
• The single-particle band problem is solved with an extension of LMTO that includes a highly precise and efficient full-potential implementation.Source: PROGRAM SUMMARY · p. 1
• An advanced fully relativistic, non-collinear implementation based on the atomic sphere approximation supports calculations of transport and magnetic properties.Source: PROGRAM SUMMARY · p. 1

1. Introduction

• The introduction distinguishes DFT basis sets by localized versus plane-wave functions and pseudopotential versus all-electron treatment; Questaal follows Slater’s all-electron augmented-wave tradition.The Slater approach joins interstitial envelope functions smoothly and differentiably to partial waves inside augmentation spheres; Questaal’s augmentation can be converged to an exact solution for the reference potential.Source: §1 · pp. 1–2, Figure 1 · p. 2

1.2. Main Features of the Questaal Package

• Questaal’s full-potential code generalises LMTO with smooth-Hankel envelope functions, offering high accuracy and efficiency, with APWs available to reach LAPW-level convergence.The package adopts the spherical approximation for basis construction; its errors are second order in the nonspherical potential component. Its recent screening transformation enables short-ranged full-potential functions and underpins the developing Jigsaw Puzzle Orbital basis, intended to combine completeness, minimality, and short range.Source: §1 · pp. 2–3, §1.2 · p. 3

1.1. Questaal’s History

• Questaal evolved from the Stuttgart code and includes early all-electron GW work; it implements QSGW+DMFT with response functions, alongside one other such implementation.The paper describes the code’s response-function implementation as unique to the authors’ knowledge.Source: §1.1 · p. 3

1.2. Main Features of the Questaal Package

• Questaal includes full-potential and ASA DFT codes, LDA+U, GW and QSGW, DMFT-solver interfaces, layer and alloy Green’s-function methods, and empirical tight binding.The ASA codes include CPA and magnetic-exchange capabilities, while the layer Green’s-function code supports transport calculations and a nonequilibrium Keldysh technique.Source: §1.2 · p. 4

1.3. Outline of the Paper

• The paper develops traditional LMTO theory, full-potential implementation, GW and QSGW, DMFT interfaces, spin and charge susceptibilities, and software aspects.It also covers coherent-potential alloy theory, noncollinear magnetism, fully relativistic techniques, LDA+U, and combined LMTO+APW methods.Source: §1.3 · pp. 4–5

2. The Muffin-Tin Potential and the Atomic Spheres Approximation; 2.1. One-centre Expansion of Hankel Functions

• KKR augments singular Hankel envelope heads with partial waves, while one-centre Bessel expansions represent remote-site tails and converge more slowly at larger radii.The head cutoff ℓb determines Hamiltonian rank; lmf typically uses ℓa = ℓb + 1, whereas traditional augmentation usually requires ℓa = 2ℓb.Source: §2 · p. 5, §2.1 · pp. 5–6

2.2. Partial Waves in the MT Spheres; 2.4. Linearisation of Energy Dependence in the Partial Waves

• Spherical partial waves encode scattering through their energy-dependent logarithmic derivatives, whose pole structure the linear method approximates over an energy window.Linearising the energy dependence of partial waves is Andersen’s central contribution to the method.Source: §2.2 · p. 6, §2.4 · p. 7

2.5. The Traditional LMTO Method; 2.6. Energy-Dependent Muffin-Tin Orbitals

• Traditional LMTO takes the envelope kinetic energy as zero, while energy-dependent muffin-tin orbitals match continuously and differentiably at sphere boundaries without diverging at large distances.In the interstitial, the orbital’s Hankel envelope is joined to a combination of a partial wave and a Bessel function inside the sphere.Source: §2.5 · p. 7, §2.6 · p. 8

2.6.1. Tail Cancellation

• Tail cancellation removes Bessel contributions inside every sphere; nontrivial KKR solutions occur when P − S has a zero eigenvalue, giving exact muffin-tin solutions in the stated limits.The limits are L → ∞ and the interstitial kinetic energy E = εi − vMTZ.Source: §2.6.1 · p. 8

2.7. The Atomic Spheres Approximation

• The ASA uses spherical site potentials, space-filling enlarged spheres with geometry violations ignored, and κ = 0 Hankel envelopes augmented by partial waves.A flat interstitial potential and a combined-correction term account for the interstitial; full-potential methods instead use better envelope functions.Source: §2.7 · p. 9

2.7.1. Tail Cancellation in the ASA

• In the ASA, linearising 1/P turns nonlinear tail cancellation into a linear algebraic eigenvalue problem and enables an efficient energy-independent variational basis without tail cancellation.Source: §2.7.1 · p. 9

2.8. Potential and Normalisation Functions

• The potential function Pℓ relates to the partial-wave logarithmic derivative and boundary conditions, while Pℓ and Nℓ are fixed by matching continuity and differentiability at the sphere boundary.Source: §2.8 · p. 9
• Linearisation of φ or Φ has second-order errors in ε − εν, giving fourth-order errors in the variational energy estimate.Source: §2.8 · p. 10

2.8.2. Spin Orbit Coupling as a Perturbation

• The formalism extends to the Pauli Hamiltonian by treating spin-orbit coupling perturbatively and including mass-velocity and Darwin terms.Source: §2.8.2 · p. 11

2.8.3. How the ASA-Tail Cancellation Reduces to a Linear Algebraic Eigenvalue Problem

• The ASA-tail cancellation condition reduces to a Hermitian linear algebraic eigenvalue problem whose eigenvalues correspond to zeros of |eP − S|.The constructed matrix is accurate to second order in ε − εν.Source: §2.8.3 · p. 12

2.9. “Screening” Transformation to Short-ranged, Tight-binding Basis

• Screening rotates the basis without changing its Hilbert space, producing short-ranged envelopes while retaining κ near zero; E=0 is a compromise because the optimal E is roughly +0.3 Ry in close-packed systems.The transformation admixes neighboring bare envelopes; the useful optimal E lies in the middle of the occupied bands, where the basis is more accurate.Source: §2.9 · pp. 12–13

2.10.1. Redefinitions of Symbols

• Screened MTO potential functions generalise the unscreened functions by matching the envelope to the partial wave and screened Bessel function; screened definitions relabel quantities by representation.The original unsuperscripted C, ∆, and p definitions correspond to the γ representation, not the unscreened α=0 representation.Source: §2.10.1 · p. 13, §2.10.2 · p. 14

2.10.2. Potential and Normalisation Functions for Screened MTO’s

• The energy derivative of the inverse screened potential function P is independent of the screening representation.The stated forms apply when P is parameterised by eP.Source: §2.10.2 · p. 14

2.10.4. Tail Cancellation in the Tight-binding Representation

• The screened MTOs exactly solve the ASA-MT potential under tail cancellation, with the practical solution rapidly convergent as the L-cutoff increases.The exact condition is subject in practice to truncating the L summations.Source: §2.10.4 · p. 14

2.10.5. Hard Core Radius

• A hard-core radius controls screening: suitably adjusted radii make the envelope function decay with distance, while radii comparable to sphere radii cause excessive damping and ringing.Andersen’s later formulation uses a hard-core radius aℓ, with nearly short-ranged functions obtainable at aRℓ=0.7sR independent of κ and ℓ.Source: §2.10.5 · pp. 14–15

2.11. MTO’s and Second Order Green’s Function

• MTOs form a complete Hilbert space for any α, which can be chosen to meet physical requirements; transformations of gγ to gβ use Eqs. (49) and (44) without depending on P parameterisation.The paper notes that short-ranged Hamiltonians motivate an empirical universal choice of α.Source: §2.11 · p. 15, §2.11.1 · p. 15

2.12. The ASA Green’s Function, General Representation

• In the ASA, G is independent of screening representation α, providing a stringent implementation test; Questaal’s lmgf and lmpg parameterise Pα to second or third order in energy.Pα can also be computed by radial Schrödinger-equation integration, with energy-dependent structure constants yielding screened KKR.Source: §2.12 · p. 16

2.13. The ASA Hamiltonian: Linearisation of the Muffin-Tin Orbitals

• Linearising energy-dependent MTOs gives an energy-independent ASA Hamiltonian; in the γ representation, second-order P parameterisation recovers the second-order ASA Hamiltonian, while third-order G is more accurate than second-order eG.Questaal’s lmgf and lmpg support either second- or third-order parameterisation.Source: §2.13 · p. 16, §2.13.3 · p. 17, §2.13.4 · p. 17

2.13.1. Potential Parameters Cα, ∆α, and oα

• In the γ representation, the basis is orthogonal apart from neglected interstitial and small p-dependent terms, and hγ eigenvalues are correct to one order higher in εi−εν than those for hα≠γ.This follows because φ and φ̇ combine in each sphere in the exact proportion at each eigenvalue of hγ.Source: §2.13.1 · p. 17

2.14. Principal Layer Green’s Functions

• Principal-layer partitioning makes the Hamiltonian block tridiagonal, so constructing G scales linearly with layer count; lmpg’s difference-equation and sparse-matrix methods scale linearly, with the former faster for large systems.lmpg constructs the auxiliary g, which can suffice for transmission and reflection, while G is needed for charge density.Source: §2.14 · p. 18, §2.14.1 · p. 18

2.14.2. Green’s Functions for the End Regions

• End-region surface Green’s functions can be found by decimation or a periodic-potential difference-equation method, but Questaal’s standard distribution lacks a fully satisfactory all-purpose method.The difference-equation method yields eigenfunctions but faces a challenging eigenvalue problem; decimation can encounter trouble where growing and decaying solutions approach unity, often at physically interesting energies.Source: §2.14.2 · pp. 18–19

2.15. Contour Integration over Occupied States

• Occupied-state Green’s-function integrals are evaluated on an elliptical complex-energy contour, with Legendre quadrature points bunched near EF; good results typically use 12–20 points.Third-order parameterisation may introduce spurious poles, which can be avoided by using the orthogonal representation and/or fairly large contour eccentricities.Source: §2.15 · p. 19

2.16. Spin-orbit Coupling in the Green’s Function

• Spin–orbit coupling is added perturbatively through matrix elements that modify the variational energy as a perturbation of potential parameters, while the basis retains definite-spin functions calculated without SOC.The energy dependence is modified by allowing the ω parameter to become a matrix.Source: §2.16 · p. 19
• For FePt, lm and lmgf magnetocrystalline-anisotropy results agree very well, and for (Fe1−xCox)2B in the virtual crystal approximation they agree essentially perfectly.The comparisons pair second-order lmgf potential functions with the two-centre lm approximation, and third-order lmgf functions with the full three-centre lm calculation.Source: §2.16 · p. 20, Figure 7 · p. 20

2.17. Fully Relativistic LMTO-ASA

• Questaal extends LMTO-ASA to fully relativistic calculations using a simplified Kohn–Sham Dirac equation that neglects the orbital contribution to the four-component relativistic current.The formulation also neglects magnetic spin-orbit coupling between radial amplitudes with different relativistic quantum numbers; the paper justifies this omission as small in most cases.Source: §2.17 · pp. 20–21

2.18. Coherent Potential Approximation

• Questaal’s CPA calculations of saturation magnetisation in Py1−xCrx track measured values fairly well, while alloy scattering broadens majority-spin bands.The CPA majority-spin spectral function also allows scattering lifetimes to be extracted.Source: Figure 8 · p. 23, §2.18 · p. 24

2.19. Noncollinear Magnetism

• At 35% Ni, Fe-Ni alloys lie near a collinear–noncollinear transition, which contributes negatively to the Grüneisen parameter.Questaal’s ASA codes support rigid-spin noncollinear configurations, spin statics, and spin dynamics, with DFT torques obtained from off-diagonal spin-density-matrix elements.Source: §2.19 · p. 24, Figure 9 · p. 24

3. Full Potential Implementation

• Questaal’s full-potential code lmf implements DFT without shape approximations and also supplies the one-body part of quasiparticle self-consistent GW.Related tools provide GW and DMFT solver inputs, and post-process quantities using the one-body part from lmf and dynamical self-energies from GW or DMFT.Source: §3 · p. 24

3.5. One-Centre Expansion

• The Jigsaw Puzzle Orbital basis tailors envelope functions to the potential, aiming for accurate, minimal basis functions and a practical realization of compactness.The current basis is a unitary transformation of the original one; compactness is particularly important in many-body treatments.Source: §3 · p. 25

3.1. Smooth Hankel Functions

• Smooth Hankel functions convolve ordinary Hankel functions with a Gaussian, making them regular near the origin while recovering Hankel behaviour far from it.Their Fourier-transform radial part is independent of angular momentum, enabling an efficient recurrence for higher-angular-momentum functions.Source: §3.1 · p. 26

3.2. Gradients of Smooth Hankel Functions

• Gradients and position operators map smooth Hankel functions into members of the same function family, enabling their matrix elements to be evaluated within that family.Gradients and energy derivatives are needed for forces, momentum matrix elements, and integration of some matrix elements.Source: §3.2 · pp. 27–28

3.3. Two-centre Integrals of Smoothed Hankels

• Products of smooth Hankel functions centred at different sites have closed-form integrals expressed as sums of other functions evaluated at the sites’ connecting vector.The result follows from Parseval’s identity and the product structure of their Fourier transforms.Source: §3.3 · pp. 28–29

3.4. Smoothed Hankels of Positive Energy

• For positive energies, the difference between conventional and smooth Hankel functions is real, although the individual functions can be complex.The supplied passage identifies this real-difference property for the s and p-related functions h0 and h−1.Source: §3.4 · p. 29

3.5. One-Centre Expansion

• The one-centre expansion enables smooth Hankels and plane waves to be combined uniformly, while allowing basis functions to be tailored to the potential.The expansion introduces a cutoff kmax, but can construct accurate minimal basis functions, including the PMT and JPO bases.Source: §3.5 · p. 29

3.6. Three-component Augmentation

• The three-component augmentation constructs matching partial-wave components and a one-centre-expanded component from an envelope function.The partial-wave combination matches the envelope function’s value and slope; general envelope functions need not be atom-centred.Source: §3.6 · p. 30
• Three-component augmentation converges more rapidly with ℓmax than conventional augmentation, while reaching the exact result as ℓmax→∞ and supporting accurate truncation at lower ℓ.The approach shares features with PAW but does not use pseudo-partial waves or projector functions.Source: §3.6 · p. 30

3.7. Secular Matrix

• The secular matrix uses Hamiltonian and overlap elements of basis functions whose envelope parts extend through space and are augmented inside spheres, with smooth and true parts matched in value and slope at the sphere boundary.Source: §3.7 · p. 31

3.7.1. Overlap and Kinetic Energy, and Output Density

• Overlap and kinetic-energy matrix elements use three-component augmented basis functions, and frozen partial waves make these elements environment-independent, suggesting a route to transferable pseudopotentials.Diagonalizing the secular matrix yields eigenfunctions and an output density formed from occupied-state basis-function bilinears.Source: §3.7.1 · p. 31

3.7.2. Core

• The frozen overlapping core approximation is reported to be better than self-consistent core treatment; a smooth Hankel component captures core-density spill-out and enables accurate shallow-core treatment without secular-matrix inclusion.This comes at the cost of strict orthogonality between those cores and valence states.Source: §3.7.2 · p. 32

3.7.3. Local Orbitals

• Local orbitals extend the energy window of valid band structures; in GW, including low-lying cores can shift valence QP levels by ∼0.1 eV, and high-lying states can shift the ZnO gap by several tenths of an eV.Many high-lying local orbitals are required for convergence in the cited ZnO case [60].Source: §3.7.3 · pp. 32–33

3.7.4. Effective Potential

• Compensating Gaussians make the electrostatic potential exact in the interstitial, while incomplete cancellation of potential components supports rapid ℓ convergence and shallow-core treatment without secular-matrix inclusion.The shallow-core treatment has minimal accuracy loss.Source: §3.7.4 · pp. 33–34

Three Component Total Energy in DFT

• Questaal calculates both Kohn–Sham and Harris total energies, with the Harris energy usually preferred because it converges faster away from self-consistency.Their difference serves as a measure of charge-density basis convergence; a larger-than-small difference indicates that a parameter such as plane-wave cutoff may be too low.Source: Three Component Total Energy in DFT · p. 34
• Three-component augmentation is more efficient because density-term differences converge faster with the angular-momentum cutoff than either term separately.The implementation also constructs the true potential correctly everywhere.Source: Three Component Total Energy in DFT · pp. 34–35

3.7.5. Forces

• The Mattheis-shift correction makes forces converge quadratically with ΔQ, and the corrected force deviation is much smaller during self-consistency iterations.The force is already reasonable at the first Mattheis-density iteration; without the correction, force convergence is only linear in the density deviation from self-consistency.Source: §3.7.5 · p. 35, Figure 12 · p. 36

3.8. LDA+U

• LDA+U uses a rotationally invariant, static mean-field treatment with double counting, while self-consistent orbital occupations can depend on their starting configuration.The simplified single-parameter form can provide orbital-dependent corrections for band gaps and defect levels, but is empirical and less accurate than QSGW.Source: §3.8 · pp. 36, 38

3.9. Relativistic Effects

• In heavy semiconductors, scalar-relativistic partial waves increasingly misestimate spin-orbit effects; a Dirac p1/2 local orbital reduces the CdTe error and lowers the PbTe gap by nearly 0.2 eV.For PbTe, the scalar-relativistic discrepancy exceeds 0.1 eV; QSGW orders the inverted L-point states correctly but somewhat overestimates the gap.Source: §3.9 · p. 39, Table 4 · p. 39

3.10. The PMT Method

• The PMT basis unifies different envelope functions, allowing plane waves and compact LMTO-like functions to be mixed without changing the total-energy, force, or density-assembly expressions.The generalized one-centre expansion supports this combination while addressing LMTO basis-completeness limitations.Source: §3.10 · p. 39

Advantages of PMT.

• Questaal’s atom-centred basis is much more compact than LAPW, but its standard moderately large basis misses about 5 mRy/atom in total energy; adding Hankel functions improves convergence.Even a small O sp basis improves LAPW convergence, and convergence continues to improve until the standard minimal Hankel basis is reached.Source: Advantages of PMT · pp. 39–40
• PMTs gain accuracy rapidly as the basis grows, with little difference between adding plane waves or more high-lying orbitals.Figure 14 compares null, hyper, small, and large basis sets as a function of the total basis size.Source: Advantages of PMT · p. 40, Figure 14 · p. 40

Drawbacks of PMT.

• PMTs can be over-complete when high-lying orbitals and plane waves are both substantial, and they complicate QSGW self-energy interpolation.The workarounds can be inefficient or finicky; PMTs are most useful as a slight augmentation of an LAPW or MTO basis, including low-cutoff plane waves for open molecular solids such as N2.Source: Drawbacks of PMT · p. 40

3.11. Floating Orbitals

• Floating orbitals improve QSGW self-energies in open systems, correcting silicon’s fundamental gap by about 0.1 eV, though their LDA effect is usually small.Questaal locates empty sites automatically; the authors expect the JPO basis to obviate the need for floating orbitals in such systems.Source: §3.11 · pp. 40–41

3.12. Screened, Short-ranged Orbitals in the Full-Potential Framework

• Screening preserves the Hilbert space while making full-potential LMTO functions short ranged, atom-centred, and defined by L character.The resulting functions can serve as projectors or minimal-Hamiltonian bases and support correlated-subspace selection; JPOs were still in development and were intended to improve compactness and kinetic-energy continuity.Source: §3.12 · p. 41
• For NiO, the screened self-energy gives a reasonably converged bandgap with first- and second-neighbour terms, and screened self-energy interpolation converges faster with k-mesh size.The two screening approaches give the same result at infinite range cutoff; screening also mitigates high-energy interpolation difficulties.Source: §3.12 · pp. 41–42, Figure 16 · p. 42
• Screened orbitals enable real-space matrix assembly and currently support calculations on a few hundred atoms per node, while QSGW interpolation can reduce the required q-point sampling.The screened basis avoids an energy threshold and loss of off-diagonal self-energy information during interpolation; the authors report equivalent results using half or, in some cases, a quarter of the original Brillouin-zone sampling in each direction.Source: §3.12 · p. 42

3.13.1. Choice of Basis Parameters

• Automatic basis setup uses free-atom fits by default or potential-based smoothing radii, while semicore states and sparse-cell cases receive additional basis functions.Semicore treatment is recommended when leaked charge exceeds 0.002e or the eigenvalue is above about −2 Ry; reference energies E1 = −0.5 and E2 = −2.0 Ry yield accurate bases for most materials.Source: §3.13.1 · p. 43, Semi-core states · p. 43, Molecular cases · p. 43

3.13.2. Comparison with LAPW Results

• Questaal’s lmf averaged 0.62 meV/atom in Delta values against WIEN2K, with agreement typical or better than other all-electron methods and especially strong for transition elements.The automatic setup was reliable; first-period and group VII molecular cases used an additional LAPW, while Au and group VI could require plane waves to converge V0 with respect to the basis.Source: §3.13.2 · pp. 43–44, Figure 17 · p. 44

4. GW and QSGW

• All-electron GW requires an auxiliary product basis spanning wave-function products, and Questaal extended pair-product evaluation of four-centre integrals from ASA to a full-potential scheme.Unlike plane waves, whose products remain plane waves, an all-electron basis cannot reuse the one-particle basis directly.Source: §4 · p. 44

4.1. Need for Self-consistency

• GW results depend strongly on the reference H0, while full scGW can worsen dielectric functions and quasiparticle levels in semiconductors and produce poor transition-metal spectra.LDA-based GW’s diagonal-only approximation can also preserve a wrong band ordering and generate unphysical crossings; QSGW’s broad reliability is not absolute, and its RPA total energies do not reach quantum-chemical accuracy.Source: §4 · p. 44, §4.1 · pp. 45, 47
• Questaal’s QSGW iterates a static Hermitian approximation to Σ to self-consistency, choosing H0 within GW rather than using full scGW’s energy-dependent, non-Hermitian self-energy.The method’s nonlocal potential includes off-diagonal components; quasiparticle self-consistency is reported to perform better than full self-consistency, though multiple competing stationary states can leave no unique optimal reference.Source: §4.1 · pp. 45, 47
• In ideal P¯3m1 TiSe2, QSGW gives a metal whereas LDA-based GW gives a small positive gap; off-diagonal Σ-induced eigenfunction and density changes explain the reversal.A density-response correction to LDA-based GW closely resembles full QSGW, including a negative gap; the passage states that the P¯3c1 phase is insulating only because of lattice displacements.Source: Figure 18 · p. 45, §4.1 · p. 46
• For MnAs and FeTe, LDA-based GW worsens the electronic structure relative to LDA, while QSGW increases MnAs exchange splitting and avoids FeTe’s unphysical flat band at EF.In MnAs, QSGW places majority and minority 3d bands at about (−2.5, 1.5) eV, versus about (−1.5, 1) eV in LDA.Source: §4.1 · p. 46, Figure 19 · p. 46

4.2. Questaal’s QSGW Implementation

• Questaal’s QSGW implementation separates lmf, lmfgwd and lmgw across the self-consistency cycle and interpolates Σ0−Vxc to obtain eigenfunctions and eigenvalues at any k.It uses a high-energy diagonal-only approximation for Σ0, typically above Ecut ∼2 Ry, and the lmf and lmgw k meshes can be independent.Source: §4.2 · p. 48, Figure 21 · p. 48

4.3. Successes of QSGW

• QSGW reliably improves many properties over DFT and LDA-based GW, but usually overestimates fundamental gaps and conduction-band effective masses.It also narrows localized d-band widths, widens widths in wide-gap semiconductors and graphene, and predicts measured nonparabolicity parameters within experimental resolution where available.Source: §4.3 · pp. 48–49
• QSGW systematically overestimates E0 and E1, while E0−E2 averages about 0.1 eV below experiment; truncating high-energy off-diagonal Σ0 contributes to this error.E2 measurements are less certain, and the cited truncation error may be revisited using a short-range basis that permits Σ0 interpolation without a finite Ecut.Source: §4.3 · p. 49, Figure 24 · p. 50
• In zincblende semiconductors, Σ-scaled QSGW almost eliminates bandgap overestimation and effective-mass underestimation, though small-gap InAs and InSb retain large relative gap errors.The hybrid scheme also improves semiconductor bandgaps and predicts Dresselhaus splitting in sp semiconductors, but is empirical and does not correct optical blue shifts or excitonic effects.Source: §4.3 · p. 50, Figure 25 · p. 50
• Without ladder diagrams, the low-frequency dielectric constant is about 80% of experiment; ladders improve charge susceptibilities and enable QSGW to describe VO2 optical conductivity.Ladder contributions account for most charge-susceptibility error, and are important to the cited VO2 result.Source: §4.3 · pp. 50–51, Figure 26 · p. 50
• For GaAs, QSGW describes Γ-L and Γ-X dispersions and valence bands well, while underestimating core-state binding and slightly overestimating conduction bands and effective mass.The nonparabolicity parameter varies by about 30% with direction, and the reported gap, effective mass, and nonparabolicity align well with experiment after the stated correction; a low-order vertex corrects much of the shallow-core discrepancy.Source: §4.3 · p. 51, Table 6 · p. 51
• With ladders and a low-order vertex correction for narrow semicore states, QSGW describes a wide range of materials uniformly, but strong spin fluctuations and rare-earth multiplet splittings remain limitations.The paper identifies occupied–unoccupied 4f splitting as too large and multiplet splittings as beyond QSGW’s scope.Source: §4.3 · p. 51

4.4. Limitations of QSGW

• In La2CuO4, nonmagnetic QSGW predicts a metal while antiferromagnetic QSGW gives a 3.5 eV gap, versus an experimental gap of about 2 eV; strong-spin-fluctuation systems expose GW’s limits.QSGW also overestimates the optical gap in La2CuO4 and CoO by about 2 eV, while ladders improve these antiferromagnetic insulators but leave discrepancies, especially in La2CuO4.Source: §4.4 · pp. 51–52, Figure 27 · p. 51
• For FeSe, nonmagnetic QSGW improves on LDA but remains farther from ARPES than elemental Fe; a low-moment SQS shifts levels closer while leaving significant discrepancies.The SQS contains three spin-up and three spin-down Fe atoms, with site moments of 0.2±0.15µB; its high-moment solution is far from experiment.Source: §4.4 · p. 52, Table 7 · p. 52
• The proposed division is to use ladders in W for reliable charge-channel susceptibilities and DMFT local spin vertices for the dominant missing spin-channel self-energy diagrams.Ladders are especially useful when included in the self-consistency cycle to construct QSGW H0.Source: §4.4 · p. 53

5. DFT+DMFT and QSGW +DMFT

• DMFT interpolates between Bloch-wave and atomic limits through a self-consistent local impurity, while its dynamic self-energy enables spectral-weight transfer beyond rigid-band descriptions.Rigid-band techniques can capture quasiparticle suppression at the Fermi energy but lack the dynamic self-energy needed to conserve spectral weight.Source: §5 · p. 53
• Questaal DMFT partitions the Hamiltonian into impurity and bath, solves the impurity, and enforces self-consistency; DFT or QSGW supplies the bath, with CT-QMC solver interfaces and.The impurity feeds back on the bath, preserving the symmetries of the entire system.Source: §5 · p. 53
• Questaal defines the correlated subspace using localized, atom-centred partial waves of a chosen angular momentum in augmentation spheres, limiting weakly correlated states in the subspace.Their localization and weaker energy dependence support an effective interaction over a wide energy window with reduced frequency dependence.Source: §5 · p. 53

5.1. Questaal’s Implementation of DMFT

• The local Hamiltonian requires a local Green’s function, screened U, and double counting; Questaal estimates U and J by cRPA, but practical calculations use static U and FLL.A fully consistent theory requires frequency-dependent U(ω), which available solvers support only in the density-density channel; excluding states over a large energy window can make U nearly energy-independent.Source: §5 · p. 54, §5.1 · p. 54
• At x=0.12, QSGW+DMFT recovers metallic La2−xSrxCuO4 with upper and lower Hubbard bands and a low-energy incoherent quasiparticle peak at EF.The calculation uses a virtual crystal approximation for x; the spectral-weight shift to low energies under doping explains the features relative to the insulating x=0 case.Source: Figure 30 · p. 54, §5.1 · p. 55
• For x=0, paramagnetic DMFT adds local spin fluctuations missing from QSGW and opens a Mott gap as Cu 3dx2−y2 splits off.The correlated subspace comprises five Cu-3d partial waves in augmentation spheres, treated with a single-site CT-QMC solver.Source: §5.1 · p. 54
• Single-site DMFT cannot capture nonlocal spin or nematic fluctuations implicated in cuprate Fermi arcs and FeSe electron-pocket suppression; cluster DMFT or Dual Fermions are alternatives.The relevant longer-ranged spin-fluctuation and nematic-fluctuation diagrams are absent from the single-site theory.Source: §5.1 · p. 55

6. Susceptibilities

• Susceptibility is a tensor response that becomes a function of one frequency and one Brillouin-zone wavevector for a time-invariant periodic reference system.In the static limit, susceptibility is the second derivative of total energy with respect to the corresponding perturbation; spin-charge coupling can make the full response a 4×4 matrix.Source: §6 · p. 55

6.1. Spin Susceptibility

• For collinear magnets, transverse and longitudinal spin responses decouple, and rotating to the x±iy basis diagonalizes the susceptibility into χ+−, χ−+, and χzz.Beyond the noninteracting response, the magnetic susceptibility requires a vertex such as the T matrix [141]; Questaal does not yet implement a perturbative T-matrix spin-susceptibility approach.Source: §6.1 · pp. 55–56, §6 · p. 55

6.1.1. Magnetic Exchange Interactions in the ASA

• The ASA Lichtenstein exchange formula is valid only in the q→0 long-wave limit and begins to deviate around k=0.25·2π/a in elemental 3d magnets.The ASA implementation derives exchange parameters using the auxiliary Green’s function and magnetic-force theorem [169].Source: §6.1.1 · pp. 56–57

6.1.2. Comparing QSGW and LDA Spin Response Functions

• QSGW predicts positive spin-wave frequencies and a Curie temperature of about 600 K for Zb-MnAs, unlike LDA’s negative frequencies throughout the Brillouin zone.The authors attribute the contrast to LDA underestimating the spin-up/spin-down 3d exchange splitting by about 1 eV, producing long-range antiferromagnetic interactions. At 25% Mn composition J remains positive but nearest-neighbour J is small; at 50% J is negative on average, suggesting a spin glass.Source: Figure 31 · p. 57, §6.1.2 · pp. 57–58

6.2. Optical Response Functions in Questaal

• Questaal’s optical tools obtain the independent-particle dielectric response from velocity matrix elements and support band-pair decomposition, core-level XAS/XES, and RIXS calculations.The approximate nonlocal velocity correction works well in weakly correlated systems but breaks down in strongly correlated NiO; the ASA second-harmonic-generation implementation has not been updated for full-potential lmf.Source: §6.2 · pp. 58–59
• Questaal’s QSGW-based Bethe-Salpeter optical spectra agree well with experiment for correlated NiO and monoclinic VO2, while the description breaks down for La2CuO4 when QSGW misses its band structure.Adding ladder diagrams to W within QSGW self-consistency is reported to improve the dielectric response even for La2CuO4; this capability is included in the standard distribution.Source: §6.2 · p. 59

6.3. Response Functions within DMFT

• DMFT-based Bethe-Salpeter calculations find dominant antiferromagnetic spin fluctuations at q=(π, π), a spin gap below a certain temperature at x=0.12, and a dx2−y2 superconducting gap.The calculations use a local vertex from CTQMC-sampled two-particle Green’s functions and non-local polarization bubbles; the spin gap is consistent with experiment.Source: §6.3 · p. 60, Figure 32 · p. 60

7. Towards a High-Fidelity Solution of the Many-Body Problem

• QSGW is the preferred non-interacting starting point and the hierarchy’s lowest level, minimally dependent on DFT and adequate for many purposes.The authors regard this as important for a consistently reliable theory.Source: §7 · p. 60
• Long-range charge fluctuations can be treated with low-order perturbation theory, whereas strong single-site spin interactions require nonperturbative treatment and weaker nonlocal contributions can be perturbative.Plasmons are typically about 5 eV and magnons about 100 meV or less; their weak interaction supports the Bose Factor Ansatz, effective in a DMFT framework.Source: §7 · p. 60
• The perturbative nonlocal route adds low-order diagrams to W, including charge-channel ladders, to improve G; electron-phonon terms are in progress and spin-fluctuation diagrams remain possible.The electron-phonon project would add a bosonic contribution to W [187]; adding low-order spin-fluctuation diagrams such as the T matrix has not yet been accomplished in Questaal.Source: §7 · p. 61
• The nonperturbative local route constructs nonlocal susceptibilities from a local interaction, then uses them to add diagrammatic contributions to G, such as Dual Fermions.The authors report remarkably good results in studied cases, including Ref..Source: §7 · p. 61
• Diagrammatic Monte Carlo could unify the routes, but realistic Hamiltonians make it prohibitively costly in time and memory; a JPO basis may enable its use for nonlocal diagrams beyond RPA.The authors expect an optimal approach to resemble the hierarchy while cross-linking its two paths, with DiagMC used for nonlocal diagrams beyond RPA.Source: §7 · p. 61

8. Software Aspects

• The code combines procedural Fortran 77 with newer Fortran features and emphasizes documentation, limited shared mutable state, argument-based data flow, performance libraries, and consistent coding practices.Development also uses version tracking and continuous integration with regression, coverage, and performance testing.Source: §8 · p. 61

8.1. Release Policy

• New changes must pass coding standards and regression tests, with at least 90% code coverage; releases progress through feature, development, and release branches.The test suite had approximately 350 tests at the time of writing; releases use vmajor.minor.patchlevel tags, and users are encouraged to update after minor-version changes.Source: §8.1 · p. 61

8.2. Parallelisation

• Questaal supports MPI and multithreading, with parallel work divided differently across its full-potential, empirical tight-binding, and GW programs.The full-potential code uses MPI mainly for Brillouin-zone sampling, while the empirical tight-binding program can distribute q-points and process groups more flexibly; GW supports multithreading and has CUDA-ported code.Source: §8.2 · pp. 61–62

8.3. Usability

• Command-line use includes input documentation and defaults, broad command-line overrides of input files, and access to online documentation, tutorials, and example scripts.Figure 34 illustrates changing k-sampling and basis options through the command line.Source: §8.3 · p. 62, Figure 34 · p. 63

8.4. Selected Publications from Questaal

• Selected publications document Questaal capabilities spanning DFT, transport, and QSGW, including self-consistent QSGW, nonequilibrium electron transport, and QSGW+DMFT+BSE susceptibilities.The section also points to publications on adiabatic spin dynamics and genetic-algorithm prediction of critical temperatures.Source: §8.4 · pp. 62–63

8.5. Distribution and Licensing

• The code is distributed under GNU General Public License version 3, and contributions under a compatible license are welcome.Historically it was distributed as source tarballs, with registered users also able to clone the online repository.Source: §8.5 · p. 63

9. Conclusions

• The paper brings classic expressions and previously unpublished developments together to explain how Questaal’s methods connect.It also summarizes the evolution of a code descended from an early all-electron method developed in Stuttgart.Source: §9 · p. 64
• The authors assess Questaal’s strengths and limits across varying levels of theory.Source: §9 · p. 64
• Questaal is presented as a promising path to efficient, high-fidelity electronic-structure calculations spanning diverse properties, materials, and levels of approximation.Source: §9 · p. 64

Appendices

• A subscript R marks site dependence, R in a function subscript makes r relative to site R, and an uppercase angular-momentum symbol such as L covers both ℓ and m.The augmentation radius is denoted s, with s_R indicating a site-dependent radius.Source: Appendices · p. 64

A. Definition of Real and Spherical Harmonics

• Questaal uses the same spherical-harmonic definitions as Jackson.Source: A. Definition of Real and Spherical… · p. 64
• The −m and +m functions are related by the relations referenced to Jackson Eqs. (3.51) and (3.53).Source: A. Definition of Real and Spherical… · p. 64
• Questaal mostly uses real harmonics Ylm, related to spherical harmonics as real polynomials in x, y, and z.Source: A. Definition of Real and Spherical… · pp. 64–65
• Products of spherical-harmonic polynomials expand as linear combinations of these functions; Gaunt coefficients vanish unless k + ℓ−m is even, leaving a polynomial in x, y, and z.Source: A. Definition of Real and Spherical… · p. 65
• For p = 1, 2, 3, the stated radial forms are energy-zero Hankel and Bessel functions; C(−) is the transpose of C(+), distinguishing ℓ+1 from ℓ−1 mappings.Source: A. Definition of Real and Spherical… · p. 65

B. Definition of Hankel and Bessel Functions

• Spherical parts use subscript ℓ, while solid functions use subscript L to denote both ℓ and m parts.This notation distinguishes radial parts of spherical functions from their solid versions.Source: B. Definition of Hankel and Bessel… · p. 65
• Questaal generally follows Methfessel’s definitions for Hankel and Bessel functions.The section distinguishes these standard definitions from Andersen’s alternative convention, marked with a caret.Source: B. Definition of Hankel and Bessel… · p. 65
• For real energy, H_L and J_L are real, and the structure constants obey S_KL(R) = S*_LK(−R).The functions can be generated from the ℓ=0 functions using Y_L(−∇); this is useful in derivations and for expressing gradients of ∇H_L as linear combinations of other H_L.Source: B. Definition of Hankel and Bessel… · p. 65
• Andersen’s caret-marked alternative definitions apply for E ≤ 0 with κ² = −E; the length scale w is arbitrary and usually averages MT radii.The alternative convention is noted as historically motivated and convenient for the ASA when E = 0.Source: B. Definition of Hankel and Bessel… · p. 65

C. Calculations in Real and Reciprocal Space

• Questaal codes work in reciprocal space except that lmpg uses real space for the principal-layer direction, and the k mesh is uniform with user-selectable centering on or symmetric straddling of k=0.The formalism, except for APWs, can apply in either real or reciprocal space.Source: C. Calculations in Real and Reciprocal… · p. 66
• Bloch sums are formed directly from short-ranged screened MTOs, while the standard lmf basis requires an Ewald technique because it is not sufficiently short ranged.Source: C. Calculations in Real and Reciprocal… · p. 66
• The lattice-vector sum is periodic in R and can be written as a Fourier series in reciprocal-lattice vectors G, yielding the Bloch-summed HL form.Source: C. Calculations in Real and Reciprocal… · p. 66

D. Matrix elements of Smooth Hankel Envelope Functions

• Parseval’s identity moves the matrix-element calculation to reciprocal space, where products of smooth Hankel functions reduce to sums from the same class.The Fourier transform of that class is known, enabling analytic expressions.Source: D. Matrix elements of Smooth Hankel… · p. 66
• Because ∇2HkL = Hk+1L, overlap and Laplacian elements belong to the same class; the overlap uses combined smoothing radius γ1 + γ2 and connecting vector R1 − R2.The derivation expands angular and radial factors before identifying the product as a Fourier transform of a linear combination of HkL functions.Source: D. Matrix elements of Smooth Hankel… · pp. 66–67
• For equal energies ε1 = ε2, the radial factor falls outside the HkL family, so the derivation introduces a new function for this limit.The special case is treated separately from the general overlap expression.Source: D. Matrix elements of Smooth Hankel… · p. 67
• The WkL functions are related to energy derivatives of HkL, and known real-space forms of H0L and its derivative provide a practical route to W0L.Higher-k WkL functions are generated from radial functions.Source: D. Matrix elements of Smooth Hankel… · p. 67
• The resulting integral also covers matrix elements between two generalized Gaussians through Hk+1L(ε=0) = −4πGkL.Its parameters k′, γ, and R are those defined for the general expression.Source: D. Matrix elements of Smooth Hankel… · p. 67

• Fourier-transform rules, together with Eqs. (116), (118), and (A.11), are used to derive the reciprocal-space action of the gradient and position operators.Source: §D.1 · p. 67
• For fixed L and p, M_L;p has two nonzero elements; the more complicated position-operator derivation is restricted to k=0.Source: §D.1 · p. 68

D.2. Matrix elements of position and gradient operators

• Gradient-operator matrix elements are linear combinations of the H_kL functions evaluated at the vector connecting atom centres.The position operator similarly involves overlap matrix elements between H_kL and ˙H_kL.Source: §D.2 · p. 68
• Position-operator integrals involving overlaps between H_kL and ˙H_kL can be evaluated in reciprocal space using Parseval’s theorem.Source: §D.2 · p. 68
• The radial part of the product involving b˙H2 is expanded and substituted into the square brackets in Eq. (D.19).Source: §D.2 · p. 68
• The derivation combines Eqs. D.25, D.26, D.21 and D.28, with Wk′M and Hk′M evaluated at the same k′, γ and R arguments as in Eq. (D.21).Source: §D.2 · p. 68