  We present an automatic, high-order accurate, and adaptive Brillouin zone
integration algorithm for the calculation of the optical conductivity with a
non-zero but small broadening factor $\eta$, focusing on the case in which a
Hamiltonian in a downfolded model can be evaluated efficiently using Wannier
interpolation. The algorithm uses iterated adaptive integration to exploit the
localization of the transport distribution near energy and energy-difference
iso-surfaces, yielding polylogarithmic computational complexity with respect to
$\eta$. To demonstrate the method, we compute the AC optical conductivity of a
three-band tight-binding model, and are able to resolve the Drude and interband
peaks with broadening in the sub-meV regime to several digits of accuracy. Our
algorithm automates convergence testing to a user-specified error tolerance,
providing an important tool in black-box first-principles calculations of
electrical transport phenomena and other response functions.
