Representative Random Sampling of Chemical Space

TL;DR

This paper introduces a method for unbiased, representative random sampling of chemical space, enabling estimation of the total number of molecules and assessing the representativeness of existing chemical databases without enumerating all molecules.

Contribution

The work presents a novel sampling technique for chemical space that is unbiased, efficient, and applicable to graph-representable molecules, along with a method to estimate total molecule counts.

Findings

The method efficiently samples molecules of up to 30 atoms.

It enables estimation of the total number of molecules in chemical space.

It assesses the representativeness of current chemical databases.

Abstract

The overwhelming majority of molecules remains unexplored. This is mostly due to the sheer number of them, which prohibits any enumeration of chemical space, the set of all such molecules. In practice, only subsets of chemical space are considered, but those subsets exhibit substantial bias, prohibiting data-driven characterization of chemical space itself. In this work, we provide a method produce unbiased representative random samples of chemical space without enumeration of constituent molecules and to estimate the number of molecules in any custom chemical space. The approach is applicable to molecules which can be represented as graph and runs efficiently even for molecules of 30 atoms. We use it to estimate the representativeness of current databases with respect to their underlying chemical space and to establish a necessary criterion for a lower bound of database sizes to be representative of an underlying chemical space.

Tables

Table 1. Table 1: Chemical spaces used in this work. Element labels are illustrative only, since any labeling of equivalently bonding elements yields the same protomolecule set.

Space	Valence	Multiplicity	Example
A	1	5	F, H, Cl, Br, I
	2	2	O, S
	3	2	N, P
	4	3	C, S, Si
	5	2	N, P
	6	1	S
B	1	5	F , H, Cl, Br, I
	2	2	O
	3	2	N, P
	4	1	C
	5	1	P

Table 2. Table 2: System size coverage of the data in this work for each chemical space. For the domain labeled exact , the correct weights are known from enumeration, so the sampling is guaranteed to be uniform. For base data, the average path length has been found heuristically from explicit sampling of all labeled degree sequences. For pure domains, this has been done for all unlabeled degree sequences only. Above the specified domains, asymptotic scaling relations are employed. base is a strict superset of pure data.

Space	Exact	Base	Pure
A	3-10 atoms	10-15 atoms	10-20 atoms
B	3-10 atoms	10-22 atoms	10-31 atoms

Table 3. Table 3: Number of molecules in the original databases (DB size) and corresponding estimated chemical space sizes, along with Kolmogorov–Smirnov (KS) and Kullback–Leibler (KL) divergence scores.

Dataset	DB size	Chemical space size	KS	KL
ANI-1	$5.75\times {10}^4$	$3.663\times {10}^6$	0.260	1.291
QM9	$1.31\times {10}^5$	$5.897\times {10}^7$	0.485	5.141
GDB-13	$9.77\times {10}^8$	$1.100\times {10}^{15}$	0.436	16.389

Equations

$$O(N_a)\equiv \{q\in p(N_a,V)|M_V(q)\text{ even}\wedge S(q)\}$$

$$C(N_a)$$

$$c(x)$$

$$l_G\sim \log |U(d)|$$

$$N_P(d)=\prod _v\prod _i\left(\genfrac{}{}{0ex}{}{{\sum }_{j\ge i}{c}_j}{c_i}\right)$$

$$l_G(d)=\left(1+{\left[\sum _i{d}_i\right]}^{-1}\log N_P(d)\right)l_G(d_U)$$

$$\log |U(d)|\approx 1.220l_G-0.7295$$

$$G=$$

$$\displaystyle\exp\Bigg{(}\left(y_{1}-\frac{1}{2}\right)\frac{M_{2}}{M}+\left(x_{2}-\frac{1}{2}\right)\frac{M_{2}^{2}}{2M^{2}}+\frac{M_{2}^{4}}{4M^{5}}$$

$$\displaystyle-\frac{M_{2}^{2}M_{3}}{2M^{4}}+\left(x_{3}-x_{2}+\frac{1}{3}\right)\frac{M_{3}^{2}}{2M^{3}}+\mathcal{O}(k_{\text{max}}^{3}/M)\Bigg{)}$$

$$M_r=$$

$$l_G\approx 0.7561\log (G)-14.40$$

$$\text{KL}=\sum _{x\in \mathcal{X}}P(x)\log \frac{P(x)}{Q(x)}$$

Full text

I Introduction

Chemical space is the set of all stable compounds from which we can choose molecules of favorable properties and therefore the search space for compound design[1, 2, 3, 4]. For practical reasons, we typically consider only a subset of that space, even when statements about general structure-property relationships are made. If this subset is not representative of chemical space, it introduces a bias which propagates to the conclusions drawn from the subset analysis. Although quantum chemistry datasets[5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19] grow larger and larger[20, 21], any self-similarity of the data points they cover does not necessarily reduce the bias of the data distribution[22, 23]. These self-similarities are also introduced by the whole tool chain that is used in generating the databases[24, 25]. Therefore, it is highly desirable to assess chemical space in an provably unbiased manner

Unbiased exploration would improve our intuition on chemical processes and modelling. Textbook concepts such as electron donating/withdrawing groups, the Hammond postulate[26] or Hückel’s rule[27] are derived from experimental data (which by the nature of the intricate experimental process is self-similar) or theoretical models. As such this might constitute overfitting towards the tiny section of chemical space we have observed so far and may break down systematically for other regions[28, 29, 30, 31], failing to generalize to an appreciable number of systems[32].

Considering main goals of modern machine learning methods, that is efficiency, accuracy, scalability, and transferability (EAST)[33], unbiased reference data connects efficiency and transferability. The smaller the required amount of training data, i.e. the higher the data efficiency, and the more evenly balanced this data is over the target domain, and the higher the transferability which can be reasonably expected.

Exploring the vastness of chemical space is relevant for a wide range of applications, such as drug design[34, 35] or catalysts[36, 37, 38], novel materials [39, 40] and molecular properties prediction [41, 42]. Given the extent of exponentially increasing chemical space, its vast majority is unexplored and probably will remain unexplored, which makes it even more important to traverse the options most efficiently. Researchers have addressed this substantial challenge by forward design based on predicting individual molecular properties with a wide range of methods[43, 44, 45, 46, 47] and inverse design using generative models[48, 49]. In forward predictions, a core challenge is to reliably assess the generalization error, which is typically done using hold out data sets, assuming that the underlying database is representative of the problem space, since otherwise aggregate measures do not apply to true out-of-sample data. For generative models, an open question is whether their image space captures the full target domain, or whether they are most likely to reproduce compounds similar to the training data or how to restrict to meaningful output spaces[50, 51]. While large training data sets were a necessary requirement for machine learning of molecular properties, most early data sets make assumptions simplifying their generation and include artificial biases towards compounds which have better support from cheminformatics toolchains[11] or quantum chemistry software[52]. Those restrictions however skip interesting systems which would commonly be assumed to be unstable, such as cubane [53], hexamethylbenzene dication [54], argon fluorohydride [55], hexanitrogen [56], C48 (cyclo[48]carbon) [57]. This impacts or hinders the discovery of of the so-called unknown unknowns, i.e. new effects or physical phenomena.

One approach of systematic exploration of chemical space is based on full enumeration. Over decades, many tools have been developed such as MOLGEN[58], OMG[59], PMG[60], surge[61], Maygen[62], and enu[63]. What they have in common is that they systematically enumerate the space of all feasible molecules by constraining themselves to those bond topologies which are in line with standard valence bond rules and then systematically generate molecular graphs. The major achievement of those codes is that duplicate checks are not required as the generation itself already guarantees uniqueness. However, even those approaches are limited in molecular size since chemical space is scaling even further than the few ten atoms that can commonly be treated. To combat the exponential scaling of chemical space, most of these tools allow the definition of further restriction such as avoiding certain features of molecular graphs or further constraining bond topology or subgroups. This has been very helpful in practical application, but only marginally increases the size of the molecules that can be enumerated systematically while introducing the aforementioned biases. In the context of structure elucidation, stochastic approaches have been pursued[64] with a focus on identifying possible target compounds if molecular fragments are known, which does not allow for random sampling of a chemical space. More recently, larger chemical diversity became the goal of molecular generators[65].

Moreover, these methods cannot be parallelized easily, even though some have been developed with parallelization in mind[60, 61]. Since in the systematic enumeration the problem cannot be subdivided evenly, even approximately random uniform sampling is infeasible with those methods, unless generating a complete list of all molecules would be feasible.

The infeasibility of full enumeration renders it necessary to employ probabilistic approaches for unbiased predictions over a chemical space constrained by domain-specific physical limits such as element variety, molecular size and stoichiometries instead of software support or practical convenience. For example, local refinement in chemical space can be carried out reliably using Monte Carlo methods[66]. However global exploration poses a considerable challenge: accessing all valid molecules requires long, sequential transitions, which are computationally expensive and difficult to parallelise. Ensuring detailed balance and unbiased sampling further complicates the process, leaving full exploration through classical Markov Chain Monte Carlo (MCMC) remains impractical. An unbiased random sampler, however, could parallelize MCMC by providing a more diverse set of starting points.

In this work, we present a method that generates approximately uniform random samples spanning a chemical space. We only require that the chemical space admits a representation of a molecule as a finite molecular graph, which excludes essentially the same application domains where vanilla SMILES representations are inapplicable, such as polymers[67]. This approach enables scalable and representative exploration of chemical space, suitable for applications in design, learning, and structural analysis.

II Methods

Our method Representative Random Sampling (RRS) produces a sequence of molecules which are approximately uniformly randomly selected from a given chemical space spanned by all molecular graphs which satisfy valence bond rules. While this is not applicable to all chemical spaces it is commonly used[67].

In RRS, we consider atoms of different valences as different atom types, e.g. phosphorus of valence 3 and phosphorus of valence 5 would be different atom types. The atom types form the set $A$ of pairs $(e,v)$ with element $e$ and valence $v$. All atoms of same valence are of the same valence type, e.g. hydrogen and all halogens. In the following, $V\equiv\{a_{2}|a\in A\}$ is the ordered set of valence types. For each valence type $v$, there are potentially multiple atom types $A(v)\equiv\{a\in A|a_{2}=v\}$, counted by the valence type multiplicity $|A(v)|$. This allows us to simplify the problem considerably, since the number of possible molecules that can be formed given a set of atoms only depends on their valence type not on their element. While enumerating graphs via graph counting methods is implemented in many codes, it is not fast enough for the problem size at hand where computing or storing the list of molecular graphs is completely infeasible. We address this by first estimating the total number of molecular graphs for each sum formula within a search space and then uniformly randomly sampling from that search space by selecting a chemical formula followed by an Markov Chain Monte Carlo sampler within that chemical formula. The key difficulty here is to reliably estimate how many molecular graphs exist for any given chemical formula.

II.1 Obtaining all chemical formulae

In principle, obtaining all chemical formulae within a chemical space is an integer partition problem where a certain number of atoms is distributed over all possible elements allowed in the chemical space. Due to the sparsity structure of those partitions, a nested approach is substantially more efficient than direct partitioning. We first find the set of constitutions $C$ for molecules in a given chemical space by repeated integer partitions. A single constitution is defined as a multiset of atom types, i.e. element and valencies.

Let $p(c,B)$ be the integer partition of the ordered set $B$ into $c$ parts $\{(x_{1},\dots,x_{|B|})|\sum_{i}x_{i}=c\}$. For each possible number of atoms $N_{a}$ in the chemical space, we obtain a set of partitions

[TABLE]

where $M_{V}(q)\equiv\sum_{i}q_{i}v_{i}$ is the sum of all degrees and $S(q)$ is true of all bonds can be saturated given the valencies[61]. $M_{V}(q)$ needs to be even such that a graph is feasible. At this stage, only the valencies are taken into account but not the valence type (e.g. so far monovalent hydrogen and fluorine would not be distinguished).

We now find all possible constitutions $C(N_{a})$ by taking the Cartesian product $c$ of all partitions within each valence type, i.e. by finding all constitutions including the atom type information

[TABLE]

This two-stage approach is computationally much more efficient yet otherwise identical to directly partitioning $N_{a}$ atoms into the corresponding elements from the chemical space, since only unique cases are generated: e.g. CF4 and CH4 correspond to the same element in $O(5)$, but become different constitutions in $C(5)$.

For each constitution, there are potentially many connected loop-free multigraphs where vertices are labeled with the chemical element. These graphs all have the same degree sequence $d$ where each vertex (an atom) is also assigned an element. To exploit symmetries (such as the CF4 case above), we only consider unique protomolecules where fictitious element labels are placeholders for all elements of compatible atom types, but including the distinction from the constitution. For example, using superscripts as valencies CH2F2 and CH2Cl2 both form protomolecules of family X4Y${}^{1}_{2}$Z${}^{1}_{2}$, and therefore are distinct from CF4 and CH4 which belong to protomolecule family X4Y${}^{1}_{4}$.

II.2 Estimating molecule counts

For small number of atoms (typically up to about 10), direct systematic enumeration of all protomolecules is feasible. In this work, we use surge[61]; all resulting counts are distributed with our software packages nablachem. Beyond this realm, we obtain approximate counts using results from network topology. First, we introduce a new graph, the universe $U(d)$, which has one vertex for each protomolecule that can be obtained from a labeled degree sequence $d$. In this universe, two vertices are connected if and only if their corresponding vertices have minimal graph edit distance. We then use the results for small world networks[68] which states that the average length $l_{G}$ of the shortest path between two vertices can be related approximately to the total number of vertices in that graph:

[TABLE]

This assumes that molecular universes can be considered a small world network which numerical results will justify. The main advantage of this approach is that we only need the average distance between elements of this universe, meaning we do not need to enumerate all of the molecules that have the same labeled degree sequence, but rather a tiny subset thereof, since the average of their pairwise distances converges quickly. This is the key step enabling the reliable estimation of molecular counts in this work. Obtaining $l_{G}$ by randomly sampling protomolecules and finding their edit distance then allows to estimate the total number of protomolecules $|U(d)|$ for that labeled degree sequence $d$. This is computationally feasible until about 20 atoms.

For any pure degree sequence, i.e. where the maximum valence multiplicity is 1, the number of distinct protomolecules is minimal amongst all degree sequences which have the same valencies but with larger multiplicities. E.g. if we have two atoms of same valency and identical element label (pure degree sequence), that yields one protomolecule. If we have two atoms of same valency and different element labels (non-pure degree sequence), we can now have three different protomolecules (with the element labels one-one, one-other and other-other). If we assume that for large and complex molecules which contain many branches and asymmetries, isomorphisms are rare, we expect that almost always each atom site is unique in chemical environment. Therefore, the number of possible protomolecules $N_{P}(d)$ upon introducing more elements with the same valency can be approximated by applying combinatorial counting to groups in the labeled degree sequence $d$:

[TABLE]

where $v$ runs over all valencies in the degree sequence, and $i$ and $j$ run over all elements of that given valency $v$. For each $i$-th element, $c_{i}$ is the number of atoms of that element. Note that for unlabeled degree sequences, $N_{P}(d)=1$. The asymptotic behavior of protomolecule counts[69] suggests that the total number of protomolecules should be inversely related to the total number of bonds $M(d)$. This yields

[TABLE]

as expression to estimate the average path length $l_{G}(d)$ for an arbitrary degree sequence if the average path length $l_{G}(d_{U})$ for the corresponding pure degree sequence is known. The number of protomolecules is then obtained from $l_{G}$ as it is for smaller molecules. This is computationally feasible until about 30 atoms.

Finally, we have reached the realm of asymptotic scaling, where formal statements about the number of protomolecules exist[69]. This allows to obtain the scaling behaviour of a quantity proportional to $l_{G}$ of any unlabeled degree sequence. Using the permutation logic of the previous paragraph allows to then obtain $l_{G}$ for any labeled degree sequence. This is extremely cheap and can be applied to arbitrarily large molecules. The bottleneck then becomes the enumeration of all sum formulas from the integer partitions in the first place.

II.3 Estimating the Minimal Graph Edit Distance

The average path length of the shortest path between two vertices in the universe is the graph edit distance, meaning the minimal sets of edits that have to be done to one molecular graph in order to obtain the other. While the graph edit distance cannot be determined in polynomial time, efficient heuristics exist. In this work, we measure the graph edit distance as the minimal Wasserstein metric between the all permutations of the adjacency matrices of two molecules. For small molecules with fewer than 100 permutations, all permutations are considered, so the minimisation is exact. For larger ones, we use the minimal Wasserstein metric obtained from five different heuristics: i) random shuffling, ii) repeatedly search for a pairwise switch of atom indices which minimizes the metric from a random initial value, iii) treating the permuation as quadratic assignment problem (QAP) with the 2-opt[70] algorithm allowing all permutations, iv) as repeated QAP with permutations freezing all but one chemical element sequentially in one scan or v) as QAP with permutations freezing all but one chemical element in multiple scans and vi) employing a depth-first graph edit distance optimizer[71]. Each such heuristic has been called about 50-100 times depending on graph size.

We empirically find that for different parts of chemical space these heuristics have different efficiency but did not analyze this further as current computational cost is acceptable.

III Results

III.1 Estimating molecule counts

The core idea to estimate the size of chemical spaces from the average path length following the small world network results in eqn. 4 first needs validation. To this end, we enumerated small molecular graphs using surge[61]. Table 1 shows the chemical spaces considered and Table 2 shows the sizes of the molecules for which there is pre-computed data available. It is important to note that the element labels are arbitrary in this work: one can easily replace one element forming a single bond with another element forming a single bond without any additional computation. The pre-computed databases are limited by how many distinct elements can be used of a given valence, the valence multiplicity.

The nablachem.space[72] python package contains all this pre-computed data and automatically chooses the most accurate method available when queried. It also contains the code to extend the databases in this work to arbitrary chemical spaces.

Fig. 2A shows the correlation between the number of molecules $n$ for a given chemical formula and the average path length $l_{G}$ as determined by our heuristic. We find that there is a strong correlation between $\log(n)$ and $l_{G}$, indicating that eqn. 4 is applicable to the data at hand. This is a key results which allows to estimate the number of molecules for a given chemical formula without any enumeration by parametrizing the proportionality from eqn. 4. When fitting to the exact counts for all 285,656 unique degree sequences covering 26,076,359,902,577 molecular graphs in our database in Fig. 2A, we obtain the following approximation

[TABLE]

We find that extending the sampling steps of our heuristic does not reduce the average path length significantly. We therefore consider the deviations in the histogram to be the consequence of the extent to which the small world network approximation is valid. Since the fit is performed on logarithmized data, it is well-balanced across the widely different number of molecular graphs in each stoichiometry: Fig. 2B shows the histogram of the relative errors made in the size estimation between the real enumerated number of graphs and the estimate from eqn. 7. Our method estimates the correct order of magnitude for each stoichiometry in 73 % of the cases. Since some stoichiometries are underestimated (the region above the linear fit in Fig. 2) and some are overestimated (the region below the linear fit), the random sampling from cumulative distributions remains balanced throughout chemical space. The strong limitation arising for our method however is that there is a systematic error for any single stoichiometry which does not average out.

Since at some point even enumerating all degree sequences of interest becomes prohibitively expensive, eqn. 5 was the consequence of a simple combinatorial argument for estimating the number of unique graphs arising from non-pure degree sequences if only the size of the pure degree sequence (or, equivalently, its average path length) is known. This hypothesis is probable, since a random molecular graph is unlikely to be symmetric, and could be confirmed in the data: Fig. 2C shows the correlation between the average path length using the explicit minimum edit graph heuristic above and the results of applying eqn. 5. Note that there is no fitting involved in this step. This greatly reduces the computational complexity from estimating the size of all degree sequences to the one of estimating the size of pure degree sequences of which there are much fewer.

Going towards even larger molecules, the heuristic determination of the average path length becomes prohibitively expensive. Fortunately, asymptotic scaling relations for the number of graphs in this case are available[69] (Theorem 1.1) and are reproduced here for simplicity:

[TABLE]

where $n$ is the length of the degree sequence $\mathbf{k}$, $[\cdot]_{r}$ denotes the falling factorial and $x_{i}=1$ if $i$ is allowed as non-loop edge multiplicity (otherwise 0) and $y_{i}=1$ if $i$ is allowed as loop-edge multiplicity (otherwise 0). Here, this means all $x_{i}=1$ and all $y_{i}=0$ for $i>0$.

Since this is an asymptotic scaling relation, the prefactor is unknown. In this work, we calibrate the scaling relation to actual data by a linear fit of $l_{G}$ to $\log(G)$ on all 148,620 pure degree sequences with more than 20 atoms for which $l_{G}$ is known (see Fig. 2) and obtain

[TABLE]

this simple model performs remarkably well and remains feasible as long as listing the integer partitions of the total number of atoms is affordable. The non-pure degree sequence count then is estimated using eqn. 5 as before.

When adding $t$ Hydrogens (or other monovalent elements) to a fixed stoichiometry, eventually the number of molecular graphs has to decrease, as the double bond equivalents of the graphs approach zero. The asymptotic scaling relation does not take this effect into account: it would grow monotonously instead, even once graphs become infeasible due to the negative double bond equivalent. Comparing to the cases for which we know the exact graph count from our database, we find that the relation tends to overestimate by $t!$, so we correct by this amount. Naturally, the extreme end of the molecule count estimates are impossible to calibrate, so we recommend to limit the use of the asymptotic scaling relation once it yields a large variance for comparable stoichiometries.

III.2 Sampling molecular graphs

Once for each chemical formula the total number of molecular graphs is estimated, we can finally sample molecules for a particular chemical formula which has been selected with a probability proportional to its size in the selected chemical space. We implement this following a Markov-Chain-Monte Carlo (MCMC) approach which allows to obtain a random multigraph of a given degree sequence[73, 74]. Since the space of molecules to sample from is restricted to those of the same degree sequence, sampling the distribution converges quickly even for large molecular graphs.

In practice, there are two main limitations: i) MCMC will find each member of the automorphism group $\textrm{Aut}(G)$ with equal probability not only each unique graph $G$ and ii) random multigraphs are not necessarily connected. We deal with the former issue via rejection sampling by accepting a potential MCMC result with probability $1/|\textrm{Aut}(G)|$ which is an acceptable inefficiency, since it only becomes costly if a degree sequence would ever be dominated by symmetric graphs. The latter problem is much more severe in practice: if for a given degree sequence the total number of connected graphs is too small, the likelihood of a random graph being connected drops rapidly. While our current implementation keeps trying and just samples this region of chemical space slowly, other random graph generation methods could potentially used to improve sampling speed of those domains. Since this does not limit sampling accuracy in all other regions of chemical space, we consider numerical efficiency of this issue a topic for future work.

The sampling code is available in nablachem.space[72] as well can be run interactively on random-molecule.org.

Given that large-scale databases of quantum chemistry data are systematically generated in a high-throughput setting, it is just natural that some artifacts of the tooling involved affect the exhaustiveness of the result. For example, the widely used QM9 database [11] which is built to almost exhaustively contain molecules from C, H, N, O, F with $\leq 9$ heavy atoms, does not contain benzene, carbon dioxide, ethylene, butadiene and other common compounds, even though they would belong to the target space of the database.

With our method we can estimate whether the given database is likely to be representative of the underlying chemical space by comparing the distribution of molecules from the underlying chemical space with the distribution of the database entries. We performed this comparison for three related molecular databases: ANI-1 [75], QM9 [11] and GDB-13 [76]. To ensure the comparison, we first recreated the chemical space corresponding to each database using our method, replicating the main compositional constraints regarding elements and atom counts. The specific conditions for each database are: ANI-1 (H, C, N, O; $\leq$ 8 heavy atoms), QM9 (H, C, N, O, F; $\leq$9 heavy atoms), and GDB-13 (H, C, N, O, S, Cl; $\leq$ 13 heavy atoms). GDB-13 has additional restrictions on stoichiometries and graph motifs to bias towards specific application domains which we explicitly did not consider here, as we aim to compare to less biased chemical spaces for novelty discovery.

All molecules from the database were grouped according to their chemical formula, resulting in frequency distributions of each stoichiometry. Sorting the stoichiometries according to their size in the reference chemical space allows to calculate a well-defined cumulative count of compounds as in Fig. 2B. Since databases are meant to form a subset, they are naturally smaller in total count than the underlying space. To account for this, the cumulative counts where normalized by dividing the cumulative counts by the total number of compounds in either space or the database. This transforms the cumulative counts into cumulative distribution functions (CDF).

Since CDF are a standard way to represent distributions, a wide array to comparison metrics exist. We evaluated the Kolmogorov–Smirnov (KS) test[77] and the Kullback-Leibler Divergence (KL)[78], as they have found applications in related works[79, 80, 81]. The KS statistic is a non-parametric test of the equality of continuous, one-dimensional probability distribution. It measures the maximum distance between the CDF of two samples. Here, those are the chemical space and the database. A lower KS value indicates that the two distributions are closely aligned over the full support, which zero being the ideal value for identical samples. The $\textrm{KL}(P\parallel Q)$ divergence instead shows how one probability distribution diverges from a second distribution.

[TABLE]

This asymmetric metric captures how much information is lost when using $Q$ to approximate $P$, where $P$ represents the chemical space and $Q$ the database distribution. We found KS and KL to agree in their assessment of the similarity of the databases and their underlying chemical space, but show results for KS only. This is because the KL metric contains undefined terms if one stoichiometry $x$ from the chemical space $\mathcal{X}$ is not present in the database, since then $Q(x)=0$ in eqn. 11, which in practice requires assuming some vanishingly small value for $Q(x)$ which then inflates the metric and makes it harder to compare between different databases of different chemical spaces. KS does not suffer from this issue and due to its normalisation to the interval [0,1] is easier to interpret.

Fig. 3 shows a) how close existing databases are to their underlying chemical space (see also Tab. 3 and b) which sampling fraction is a lower bound to representatively capture the distribution of stoichiometry diversity in a database. It is not only helpful to see how well a database represents the underlying target space, but also by how much a database could be further reduced while still achieving the same representativeness compared to the full database. Smaller databases would allow for higher quality reference calculations and therefore are generally more desirable. We achieve this by random subsampling of the databases in question.

For the QM9 and GDB13 datasets, the KS values are higher than for ANI-1. This behaviour is expected given the broader chemical diversity and larger size of these datasets. Higher KS scores reflect the inclusion of valid stoichiometries that are not represented in the original datasets. The small difference between QM9 and GDB13 is the consequence of them being related: QM9 is obtained by further filtering GDB13 to those molecules which converged to a DFT minimum using a particular workflow, so naturally some of the worse score of QM9 can be explained by the selection process of which molecular graphs are stable. As discussed above, this only explains part of the QM9 score as that database is also missing practically relevant molecules. KL scores allow no direct comparison between databases of different chemical spaces as the KL divergence is particularly sensitive to mismatches in the tails of the distributions[80], which affects non-sampled low-probability stoichiometries. These outliers increase the referent distance without necessarily compromising the representativeness of the core distribution.

Interestingly, the generated chemical spaces required significantly smaller subset sizes to approximate the overall distribution with high fidelity. As a result, the generated spaces can be characterized more efficiently with fewer samples, a useful property for downstream applications such as virtual screening [82] or training data selection [22].

IV Conclusions

Whenever we make data-driven statements in an attempt to characterize global behavior of chemical space such as structure property relationships, then the fundamental question is whether those are valid for single regions in chemical space or whether this has general predictive power. Either of them would be useful, but it is important to be able to distinguish the cases. Our method potentially allows to put any such data-driven statement on a more rigorous footing by sampling molecules from a given chemical space, making sure to consider all of them without rejection and then testing whether a statement is correct on the sampled molecules. Potentially, this could either uncover a bias in existing databases or generative model output and analysis or prove the lack thereof. In both cases, this enables work towards improved understanding of chemical space as a whole, moving away from considering constituent molecules on a case by case basis. Statistical analysis of the generated molecules may help to identify additional trends and structural properties of chemical space.

One of the main advantages of machine learning is the existence of formal guarantees on how a model must behave or must improve if additional data points are provided. These formal guarantees, however, usually require a certain well-behavedness of the underlying target function, such as independent and identically distributed training data. It is not clear whether current databases fulfil that requirement since ensuring a certain distribution of data points is rarely enforced in database generation. Coverage of chemical space becomes more important with multi-level machine learning methods and with our approach, relevant coverage becomes computationally feasible even for large spaces.

The main limitations of our approach lie in a) the need for a pre-built database of average path lengths, even it is shared amongst all researchers and b) the approximative nature of the exploited correlations which renders the resulting samples only approximately uniformly distributed over chemical space. Finally, some of the resulting graphs will not admit a stable minimum energy geometry; future work may address this by introducing additional weights for stoichiometries.

Finally, our work helps to address a more subtle effect regarding machine learning databases: many of them are quite old (more than ten years[11]) but are still used for benchmark. This means that over time there is a risk that we do not obtain new models that are of lower generalization error, but rather models which, if trained on these databases, can predict the same benchmark databases well. This is a necessary consequence of using the same benchmark databases for many iterations of different models, even if each of which is designed and cross-validated following best practices. This means that in order to test generalization error, we also need test data to be a moving target as guaranteed to be unseen data[83]. Using the randomized sampling approach, this can be achieved easily because the space to sample from is virtually unlimited. In this approach, one could in regular intervals obtain small sets of unseen data that then affords a more truthful approximation to the generalization error in quantum chemistry.