How do you count molecules that don't exist yet?

1A molecule is a graph

For counting purposes, a molecule's structural formula is all that matters: which atoms are bonded to which, and with what bond order. In mathematical terms that is a graph. Atoms are vertices labeled by element, and bonds are edges labeled 1, 2 or 3.

Three conventions keep the counting well defined. They are the conventions of the GDB databases (GDB-11, -13, -17), so the numbers can be compared with them.

• Size is the number of heavy atoms, meaning atoms other than hydrogen. Hydrogens are not placed explicitly: every atom is filled up to its normal valence with H. Ethanol, CH₃CH₂OH, has 3 heavy atoms.
• Only constitutional isomers count. Two molecules differ if their atoms are connected differently. Mirror images and cis/trans forms (stereoisomers) count once.
• Molecules are neutral, closed-shell and connected. There are no ions, no radicals and no salts.

With carbon only and 2 heavy atoms there are exactly three molecules:

With carbon and oxygen there are two more: methanol (CH₃OH) and formaldehyde (CH₂=O). The O–O molecule (hydrogen peroxide) is excluded by a rule explained below, so the count is 5.

2Why naive counting fails

A first idea is to number the atoms 1…n, decide for each pair whether they are bonded, and count the valid choices. This overcounts badly, because a molecule's atoms carry no numbers. Propane is the same molecule whether its central carbon is called atom 1, 2 or 3. Numbering the atoms gives 3 versions of the propane skeleton and only 1 of the cyclopropane triangle. For larger molecules the factor approaches n! and depends on each molecule's symmetry, so it cannot simply be divided out.

The usual remedy is a canonical form: a procedure that renumbers any graph in a standard way, so that isomorphic graphs come out identical. Then "same molecule?" reduces to "same canonical form?", which a hash table answers instantly. Our canonical labeling is written from scratch in pure Python. It works in two stages:

1. Color refinement. Every atom starts with its element as its color. The program then repeatedly re-colors each atom by (its own color, the multiset of its neighbors' colors) until the colors stop changing. After this, the atoms of propane split into "end carbon" and "middle carbon", and in most molecules every atom ends up with a distinct color.
2. Individualization. If several atoms still share a color (in a symmetric molecule such as benzene, every carbon does), the program tries giving each candidate a unique color in turn and refines again. It keeps the numbering that gives the lexicographically smallest adjacency matrix. Symmetries found along the way prune equivalent branches.

3Step 1: skeletons, each exactly once

The work is split into two layers. First the program builds skeletons: connected graphs whose atoms are labeled by element but whose bonds have no order yet. It then decorates each skeleton with bond orders (step 2). One skeleton, such as a chain of three carbons, supports several molecules (propane, propene, propyne, allene).

Skeletons are grown one atom at a time. Each new atom bonds to one or more existing atoms. Grown carelessly, the same skeleton appears many times, once for each order in which its atoms could have been added. The program avoids this with canonical augmentation, a technique due to Brendan McKay (1998):

Each child therefore has exactly one accepted parent, and the skeletons form a tree in which every skeleton appears exactly once. No global list of "seen" structures is needed, so the tree can be explored in pieces, which matters for the sampling in step 3. The figure shows the start of the all-carbon tree. The 3-atom chain is the parent of the 4-atom chain, the branched star and the 4-ring. The triangle is the parent of the three skeletons that contain a triangle.

The number of carbon skeletons at each level, 1, 1, 2, 6, 21, 78, 353, 1929, … equals the number of connected graphs with every vertex of degree ≤ 4. The program reproduces this known sequence, which is one of its correctness checks.

Some rules can be checked on skeletons, which prunes the tree early. A rule qualifies if it holds for every fragment of a valid molecule:

• no atom has more neighbors than its valence (carbon at most 4, oxygen at most 2, fluorine 1);
• no forbidden element pairs are bonded (O–O, N–O, halogen–heteroatom; see the rules below).

Pruning with such hereditary rules is safe: if a skeleton breaks one, every larger skeleton grown from it breaks it too.

4Step 2: bond orders and hydrogens

For each skeleton, the program tries every combination of bond orders (1, 2 or 3 on each bond). The number of hydrogens on each atom is whatever completes its valence:

An optional --hypervalent mode also allows S with valence 4 or 6 and P with valence 5, as in sulfones and phosphates.

Resonance counts once. Benzene can be drawn with its double bonds in two alternating positions (Kekulé structures). Both drawings describe one molecule. The program therefore identifies a molecule by its skeleton plus the hydrogen count on each atom, not by the exact placement of double bonds. Both benzene drawings have one H on every carbon, so they collapse into one entry. The molecule is accepted if any of its bond-order assignments passes the rules.

5The stability rules

Valence alone allows many structures that no chemist would expect to keep in a bottle: strained rings with triple bonds, chains of oxygen atoms, molecules that rearrange immediately into something else. The program applies four groups of rules, modeled on the filters used for the GDB databases. Each group can be switched off from the command line, so you can see what it does.

Rules in action: three-carbon molecules

Valence allows 9 molecules with three carbons. The strain rules remove the four that are 3-membered rings containing multiple bonds, which leaves 5.

Rules in action: heteroatoms

Each of these structures satisfies valence, but a heteroatom filter rejects it:

6How much each rule removes

The table gives exact counts at each size, adding one rule group at a time (each row includes the rules above it). Choose an element set:

Two things stand out. The rules remove a large fraction of structures, and the fraction grows with size: for C, N, O, F at 6 atoms only about 21% of the valence-allowed molecules survive. Still, the numbers grow roughly geometrically, by a factor of 5–10 for each added atom with mixed elements. This geometric growth is the "combinatorial explosion" of chemical space.

7Step 3: estimating what you can't list

Exact enumeration in pure Python is practical up to about 7 heavy atoms (for C, N, O, F: about 1 s for 6 atoms and about 10 s for 7). The optional Rust core is about 20× faster per core and uses all cores, which buys 1–2 more atoms: 8 atoms take about 1.5 s on 10 cores. At 17 atoms there are up to about 10¹⁷ molecules, too many to list in any reasonable time. Because each skeleton appears exactly once in the tree, the program can instead estimate the size of the tree by sampling it.

Knuth's estimator: one random walk

Donald Knuth (1975) proposed a simple way to estimate the size of a tree you can't afford to traverse. Start at the root and repeatedly step to a randomly chosen child, writing down the number of children b₁, b₂, … at each step. The product

is an unbiased estimate: averaged over many random walks, it gives the true count. Multiplying by the number of molecules on the skeleton reached at each depth estimates the molecule count instead. The catch is that the estimate is heavy-tailed. Most walks wander into sparse corners of the tree, and a rare walk through a bushy region carries an enormous weight. In our tests, runs with 300 and 1000 walks to 10 atoms disagreed by a factor of 3, and the error bars were not reliable. Knuth's estimator is still available (--method knuth), but it is not the default.

Stratified sampling: many walks that share the work

The default is Pang-Chi Chen's (1992) stratified sampling, which runs many coordinated walks at once:

1. Start with every skeleton at the exact frontier (6 atoms), each with weight 1.
2. Expand all current skeletons into their children. Each child inherits its parent's weight.
3. Group the children into strata, classes of skeletons expected to have similar subtrees. Here a stratum is "how many atoms have 0, 1, 2, … free valences, and how many rings". Within each stratum, keep one representative, give it the total weight of the group, and choose which member survives at random with probability proportional to its weight. This keeps the estimate unbiased.
4. The estimate for this size is the sum over representatives of weight × (molecules on that skeleton). Then repeat from step 2 for the next size.

Each pass explores a broad cross-section of the tree, one representative of every structural type, instead of a single path. It costs more than one Knuth walk but has far lower variance. Running 8–10 independent passes gives a mean and a standard error:

The relative error grows slowly with size. At 17 atoms, the per-size estimates in the explorer's data (12 passes per element set) carry a standard error of about 4–34%, depending on the element set; the largest sets are the hardest. Totals over all sizes come out at 3.5–32%. The explorer shows a ±2 SE band.

8Beyond the computed sizes

The data behind the explorer is computed up to 17 heavy atoms. For larger sizes the explorer extrapolates, labels the result as an extrapolation, and shows an envelope rather than a single line:

• Upper edge: the growth factor per added atom stays at its value between 16 and 17 atoms.
• Lower edge: the growth factor decays gradually toward the alkane growth rate. Branched alkanes grow by a factor of about 2.95 per carbon in the limit of large size. The lower edge is never allowed to fall below the alkane count itself.
• Central value: the geometric midpoint of the two edges.

The upper edge is not a true bound. The computed growth factor per atom keeps rising with size: for pure carbon it goes from 5.4 at 8 atoms to 13 at 17, and for C, N, O, S plus halogens from 10 to 18.5. Larger molecules can close more rings, so each added atom opens more new possibilities than the one before. Holding the last factor constant therefore tends to underestimate. When the computed range was extended from 12 to 17 atoms, the GDB-like preset at 20 atoms rose from about 4×10¹⁸ to about 3×10²⁰, above the earlier envelope. Treat the large-size numbers as order-of-magnitude estimates at best. The next section compares them with published estimates.

For a direct comparison at a computed size: with the GDB-17 elements (C, N, O, S and the halogens) we count about 1.3×10¹⁷ molecules up to 17 atoms. GDB-17 lists 1.66×10¹¹, about a million times fewer, because its functional-group, ring and stability filters are far stricter than our basic rules.

9Other estimates of chemical space

Published estimates of the number of possible organic molecules span more than 50 orders of magnitude. They are not in conflict so much as answering different questions. Each one fixes a size limit (heavy atoms or molecular weight), a set of elements, and a notion of which structures count: any graph that obeys valence, molecules that pass stability filters, or molecules that can actually be made from available reagents.

What separates the numbers. Three things move an estimate by many orders of magnitude.

• Size. Counts grow faster than exponentially with the number of atoms (section 8), so the size limit dominates. The Polishchuk fit alone gives about 10²⁶ at 30 atoms and 10³³ at 36. Comparing 10⁶⁰ at 30 atoms with 10¹¹ at 17 says little about the methods behind them.
• Which structures count. Bohacek's figure is a rough combinatorial count with no stability filtering. The authors call their 10⁶³ "a very approximate estimate" and conclude only that the true number is "well in excess of 10⁶⁰". They also expect stereoisomers to raise the count, whereas GDB and this explorer count constitutional isomers only. GDB-17 removes many valid but reactive or strained structures, and its filters are much stricter than this explorer's basic rules. At the same 17 atoms and the same elements, our count is about a million times larger than GDB-17 (section 8).
• Whether anyone could make them. Ertl's estimate and make-on-demand libraries such as Enamine REAL ask a different question: how many molecules can be made with known reactions and available reagents. These numbers are much smaller than the space of stable structures, but every molecule in them is, in principle, orderable.

Where this explorer fits. Like GDB, it counts exhaustively (exactly, then by unbiased sampling) and applies stability rules, but its rules are looser, so its numbers sit between GDB-17 and the purely combinatorial estimates. The Polishchuk extrapolation is a fit to the GDB-17 counts, so it carries GDB's strict filters to larger sizes. Our growth factor per atom keeps rising up to 17 atoms, which a simple fit of that form may not capture, so extrapolations to 30 or 36 atoms remain uncertain by many orders of magnitude, ours included.

1. R. S. Bohacek, C. McMartin & W. C. Guida, "The art and practice of structure-based drug design: a molecular modeling perspective", Med. Res. Rev. 16 (1), 3–50 (1996); the estimate is footnote a to Fig. 6, p. 43. doi:10.1002/(SICI)1098-1128(199601)16:1<3::AID-MED1>3.0.CO;2-6
2. P. G. Polishchuk, T. I. Madzhidov & A. Varnek, "Estimation of the size of drug-like chemical space based on GDB-17 data", J. Comput.-Aided Mol. Des. 27 (8), 675–679 (2013). doi:10.1007/s10822-013-9672-4
3. P. Ertl, "Cheminformatics analysis of organic substituents: identification of the most common substituents, calculation of substituent properties, and automatic identification of drug-like bioisosteric groups", J. Chem. Inf. Comput. Sci. 43 (2), 374–380 (2003). doi:10.1021/ci0255782
4. L. Ruddigkeit, R. van Deursen, L. C. Blum & J.-L. Reymond, "Enumeration of 166 billion organic small molecules in the chemical universe database GDB-17", J. Chem. Inf. Model. 52 (11), 2864–2875 (2012). doi:10.1021/ci300415d
5. O. O. Grygorenko, D. S. Radchenko, I. Dziuba, A. Chuprina, K. E. Gubina & Y. S. Moroz, "Generating multibillion chemical space of readily accessible screening compounds", iScience 23 (11), 101681 (2020). doi:10.1016/j.isci.2020.101681

10How we know it works

• Alkanes. With carbon only, no rings and no multiple bonds, the program reproduces the known numbers of alkane isomers: 1, 1, 1, 2, 3, 5, 9, 18, 35, 75 for 1–10 carbons (OEIS A000602).
• Graph counts. The carbon skeleton counts match the number of connected graphs with maximum degree 4.
• Brute force. For up to 4 atoms, the tree is compared with a slow enumeration of every numbered graph, deduplicated by canonical form. They agree, which confirms that every skeleton appears exactly once.
• Hand counts. C₂ gives 3 and C,O at 2 atoms gives 5. Benzene counts once. Each rule removes exactly the structures it should.
• Samplers vs exact. Both estimators fall within 4 standard errors of the exact count for small sizes.
• GDB-11 sanity check. With C, N, O, F up to 11 atoms we estimate about 6×10⁸ molecules, about 25 times the 2.6×10⁷ in GDB-11. GDB uses much stricter functional-group and ring filters, so a larger count is expected. The same order of magnitude shows that the rules are doing roughly the right job.

11What the numbers do and don't mean

• These are candidate structures that pass simple rules, not known or synthesizable compounds. For comparison, the CAS Registry lists a few hundred million substances.
• Tautomers count separately (apart from the enol ban). The strain and Bredt rules are heuristics, not energy calculations.
• Charged species, radicals, stereoisomers and most inorganic chemistry are not counted.
• Stricter rules, such as GDB's functional-group filters or drug-likeness, would give smaller numbers. Looser ones (hypervalent S and P, N–O bonds) give larger ones. The ranking between element sets and the rate of growth are more robust than any single number.

12Try it yourself

The code is pure Python with no dependencies: github.com/Autonomous-Scientific-Agents/Chemical-Space-Explorer.

git clone https://github.com/Autonomous-Scientific-Agents/Chemical-Space-Explorer
cd Chemical-Space-Explorer
python -m chemuniverse --max-atoms 8 --elements C          # exact up to 6, sampled beyond
python -m chemuniverse --max-atoms 6 --elements C N O --no-strain     # see what strain rules do
python -m chemuniverse --max-atoms 10 --elements C N O F --no-rings --no-multiple-bonds
python -m chemuniverse --max-atoms 11 --elements C N O F --samples 10 --seed 1

The program prints the count at each size, its standard error, whether it is exact or estimated, and the running total.