  Hybrid density functional (HDF) approximations usually deliver higher
accuracy than local and semilocal approximations to the exchange-correlation
functional, but this comes with drastically increased computational cost.
Practical implementations of HDFs inevitably involve numerical approximations
-- even more so than their local and semilocal counterparts due to the
additional numerical complexity arising from treating the exact-exchange
component. This raises the question regarding the reproducibility of the HDF
results yielded by different implementations. In this work, we benchmark the
numerical precision of four independent implementations of the popular
Heyd-Scuseria-Ernzerhof (HSE) range-separated HDF on describing key materials'
properties, including both properties derived from equations of states (EOS)
and band gaps of 20 crystalline solids. We find that the energy band gaps
obtained by the four codes agree with each other rather satisfactorily.
However, for lattice constants and bulk moduli, the deviations between the
results computed by different codes are of the same order of magnitude as the
deviations between the computational and experimental results. On the one hand,
this means that the HSE functional is rather accurate for describing the
cohesive properties of simple insulating solids. On the other hand, this also
suggests that the numerical precision achieved with current major HSE
implementation is not sufficiently high to unambiguously assess the physical
accuracy of HDFs. It is found that the pseudopotential treatment of the core
electrons is a major factor that contributes to this uncertainty.
