  The work function is the key surface property that determines how much energy
is required for an electron to escape the surface of a material. This property
is crucial for thermionic energy conversion, band alignment in
heterostructures, and electron emission devices. Here, we present a
high-throughput workflow using density functional theory (DFT) to calculate the
work function and cleavage energy of 33,631 slabs (58,332 work functions) that
we created from 3,716 bulk materials, including up to ternary compounds. The
number of materials for which we calculated surface properties surpasses the
previously largest database, the Materials Project, by a factor of $\sim$27. On
the tail ends of the work function distribution we identify 34 and 56 surfaces
with an ultra-low (<2 eV) and ultra-high (>7 eV) work function, respectively.
Further, we discover that the $(100)$-Ba-O surface of BaMoO$_3$ and the
$(001)$-F surface of Ag$_2$F have record-low (1.25 eV) and record-high (9.06
eV) steady-state work functions without requiring coatings, respectively. Based
on this database we develop a physics-based approach to featurize surfaces and
use supervised machine learning to predict the work function. We find that
physical choice of features improves prediction performance far more than
choice of model. Our random forest model achieves a mean absolute test error of
0.09 eV, which is more than 6 times better than the baseline and comparable to
the accuracy of DFT. This surrogate model enables rapid predictions of the work
function ($\sim 10^5$ faster than DFT) across a vast chemical space and
facilitates the discovery of material surfaces with extreme work functions for
energy conversion, electronic applications, and contacts in 2-dimensional
devices.
