An example of the average pair distribution functions 𝑔 LiO ̅̅̅̅̅̅(𝑟) and 𝑔 LiLi ̅̅̅̅̅̅(𝑟) built for Li-O and Li-Li pairs using the ICSD data and the proxy potentials produced from them is shown in Fig. . The process, described in detail in Supporting Information, starts from selecting the structures that contain the atoms of interest in the targeted atomic environments, compiling the frequency histograms for pair separations in each structure and applying Gaussian smoothing. These histograms are then converted into radial distribution functions and combined across all reference structures, here as a simple arithmetic mean, to produce the average radial distribution function 𝑔 𝛼𝛽 ̅̅̅̅̅(𝑟). Figure shows the 𝑔 𝛼𝛽 ̅̅̅̅̅(𝑟) and its -ln𝑔 𝛼𝛽 ̅̅̅̅̅(𝑟) transform for the "raw" ICSD data without any smoothing applied in order to contrast them with the SPPs that are smooth and differentiable at any point while still closely following the raw data. In addition to the Gaussian smoothing, we used Savitzky-Golay smoothing to extract the characteristic low-frequency features from noisy 𝑔(𝑟) data typically observed for a small number of reference structures. Padding at short distances and tapering at a cut off distance 𝑟 cut = 10 Å was implemented to make sure that the SPPs produced strong repulsion for atomic overlaps and asymptotically approached zero at 𝑟 cut . The complete details and the link to the statistical data derived for all 3,388 pairs of atoms reported to co-exist in the ICSD structures are in Supplemental Material.