Supplementary information: Discovery of stable surfaces with extreme work functions by high-throughput density functional theory and machine learning

The custom algorithm to determine all unique surface terminations is illustrated in Fig. 1 of the main paper and described in more detail here. First, a list of nearest neighbors (considering only atoms above the reference atom, i.e. with a larger $c$-component, up to a cutoff radius of 10 $\AA{}$) with their distances and chemical elements are compiled for each atom in the slab. Atoms with identical nearest neighbor lists are grouped together into local environment groups (LEGs). In the next step, for each atomic layer (i.e., atoms with the same $c$-component considering a tolerance of 0.1, 0.5, or 1 $\AA{}$, depending on which yields the lowest number of unique terminations) the LEG of all atoms in the layer is grouped into a list (one list per layer of atoms). Lastly, the number of unique lists of LEGs $n_{\mathrm{env}}$ determines the number of unique terminations $n_{\mathrm{term}}$ as follows: $n_{\mathrm{term}}=n_{\mathrm{env}}$ if the slab is symmetric (i.e., a mirror or glide plane parallel to the surface, or a 2-fold rotation axis normal to the surface), or $n_{\mathrm{term}}=2n_{\mathrm{env}}$ else. This is due to the fact that the local environments are determined in the positive $c$-direction (only above the reference atom). Hence, a termination A on one side of the slab with termination B underneath might not be equivalent to the A termination on the other side if the slab is not symmetric (as illustrated in Fig. S2).

In a previous version of the database we included materials with a non-zero band gap during screening and trained machine learning models on it. However, eventually this database was not utilized in the final paper as band gaps for PBE-level DFT calculations are well-known to be unreliable and hence the work function (positioned at the center of the gap) may also suffer from the same limitation. We do however note that this error is potentially small for materials with a small band gap (less than about half an eV) with a potential partial error cancellation w.r.t. the vacuum energy level. A demonstration of that is the fact that we were able to identify low work function surfaces in non-metallic materials such alkali oxides (e.g., Cs₂O₂), alkaline vanadates (e.g., BaVO₃), and alkali-based anti-fluorite structures (e.g., Rb₂O, Rb₂S) that were already known in literature (cf. introduction in the main paper). The lowest work function after ionic relaxation was calculated for the (100)-Rb surface of PdRb₂C₂ at 0.93 eV.
Most high work function candidates in the previous version of the database exhibit a larger bandgap, as expected for flourides, carbides, nitrides and oxides. However, there are a few materials with a zero PBE bangdap in spacegroups 66 (AgO), 139 (CsO₂, RbO₂, KO₂, CaN₂), 221 (ReO₃), and 225 (PbO₂, BiO₂). The highest work function after ionic relaxation was calculated for the (110)-F surface of AlF₃ at 11.04 eV, likely not in agreement with experimental work functions due to the large band gap error in PBE.
The previous version of the database and the machine learning models trained on that database are available at https://github.com/peterschindler/WorkFunctionDatabase .