Woptic: optical conductivity with Wannier functions and adaptive k-mesh refinement

We present an algorithm for the adaptive tetrahedral integration over the Brillouin zone of crystalline materials, and apply it to compute the optical conductivity, dc conductivity, and thermopower. For these quantities, whose contributions are often localized in small portions of the Brillouin zone, adaptive integration is especially relevant. Our implementation, the woptic package, is tied into the wien2wannier framework and allows including a many-body self energy, e.g. from dynamical mean-field theory (DMFT). Wannier functions and dipole matrix elements are computed with the DFT package Wien2k and Wannier90. For illustration, we show DFT results for fcc-Al and DMFT results for the correlated metal SrVO$_3$.

💡 Research Summary

The paper introduces a novel adaptive tetrahedral mesh refinement algorithm for Brillouin‑zone (BZ) integration in three‑dimensional crystalline solids and implements it in the open‑source package woptic. Traditional BZ integration relies on a uniformly spaced k‑mesh, which is inefficient for response functions such as optical conductivity, dc conductivity, and thermopower that receive dominant contributions from limited regions of the BZ (e.g., near the Fermi surface). The authors address this inefficiency by dynamically refining tetrahedra where the estimated integration error is large, while preserving mesh regularity and shape stability.

The algorithm starts from a coarse tetrahedral mesh generated by a Kuhn triangulation of the reciprocal unit cell. Each tetrahedron belongs to one of two classes; refinement splits a tetrahedron into eight sub‑tetrahedra in a way that guarantees that the new elements remain within the same two classes, thus controlling shape distortion. Hanging nodes that appear on edges during refinement are eliminated through a mesh‑closure step, which recursively refines neighboring tetrahedra to maintain a smooth resolution gradient across the BZ. This prevents abrupt changes in mesh density that could otherwise cause numerical instability.

Error estimation is performed on a per‑tetrahedron basis. The integrand of the optical conductivity, denoted g(k), is approximated inside each tetrahedron using a low‑order Gaussian quadrature (four‑ or five‑point). The difference between the quadrature estimate and the exact integral over the tetrahedron provides an error indicator ε_T. The total error ε_tot is the sum of ε_T over all tetrahedra. The refinement loop proceeds until ε_tot falls below a user‑defined tolerance, thereby concentrating k‑points only where the integrand varies rapidly (e.g., near sharp features induced by a small imaginary part of the self‑energy).

The physical formulation of the conductivity is expressed in a Wannier‑function basis. Starting from DFT Bloch states, maximally localized Wannier functions are constructed with Wannier90, yielding a unitary transformation U(k). The Hamiltonian in Wannier space is H(k)=U†(k) E(k) U(k), where E(k) contains the DFT eigenvalues. A local self‑energy Σ(ω) obtained from a DMFT impurity solver is added to H(k) to form the interacting Green’s function G(k,ω)=