 Software
 Open Access
 Published:
Surge: a fast opensource chemical graph generator
Journal of Cheminformatics volume 14, Article number: 24 (2022)
Abstract
Chemical structure generators are used in cheminformatics to produce or enumerate virtual molecules based on a set of boundary conditions. The result can then be tested for properties of interest, such as adherence to measured data or for their suitability as drugs. The starting point can be a potentially fuzzy set of fragments or a molecular formula. In the latter case, the generator produces the set of constitutional isomers of the given input formula. Here we present the novel constitutional isomer generator surge based on the canonical generation path method. Surge uses the nauty package to compute automorphism groups of graphs. We outline the working principles of surge and present benchmarking results which show that surge is currently the fastest structure generator. Surge is available under a liberal opensource license.
Introduction
Chemical structure generators enumerate or generate molecular graphs of organic or bioorganic molecules. They are an integral part of systems for computerassisted structure elucidation (CASE) [1] and can be used to create molecular libraries for virtual screening [2, 3] or enumerate chemical spaces in general [4]. The history of chemical graph generators goes back at least to the 1960s DENDRAL project which was aimed at the CASE of organic molecules based on mass spectrometric data [5]. DENDRAL was developed for NASA’s Mariner program to search for life on Mars [5, 6]. Its structure generator used substructures as building blocks and was able to deal with overlapping substructures. In the early history of the structure generators, ASSEMBLE was another building block based structure generator [7]. In the field, there is a family of generators based on mathematical theorems such as algorithmic group theory [8] and combinatorics [9]. Besides DENDRAL, MASS [10] was also another good example for the applications of mathematical theorems in structure generation. It was a tool for the mathematical analysis of molecular structures. SMOG [11] was the successor of the MASS algorithm.
We have recently reviewed the history of chemical graph generators in detail [12].
While most structure generators work in a deterministic way, i.e. exhaustively generate structures according to given boundary conditions [13], stochastic generators were also suggested for large molecular spaces [14].
Among the currently available structure generators, such as DENDRAL, ASSEMBLE, SMOG, COCON [15] and LSD [16], MOLGEN [17] constituted the stateoftheart for decades in terms of speed, completeness and reliability.
The first version of MOLGEN was based on the strategy of DENDRAL software and developed to overcome the limitations of DENDRAL [18]. The software is based on the orderly graph generation method [19]. Although MOLGEN is the de facto gold standard in the field, it has the downside of being closedsource software. The algorithm cannot be further developed or modified by scientists based on their interests. The most efficient and fast opensource chemical graph generator was MAYGEN [20] based on the orderly generation method. However, MAYGEN is approximately 3 times slower than MOLGEN.
The state of the art of large scale structure generation was recently set by the lab of JeanLouis Reymond [21] in developing an inhouse solution for a nautybased structure generator, which enabled them to produce the numeration of 166 billion organic small molecules in the chemical universe database GDB17. To the best of our knowledge, this inhouse generator was not released as opensource or otherwise.
Thus, there is still the need for an efficient opensource chemical graph generator. In [20] we expressed the hope to “trigger a surge in the development of improved and faster” structure generators. Here we present the novel structure generator surge, based on the principle of the canonical generation path method. Surge is opensource and outperforms MOLGEN 5.0 by orders of magnitude in speed. Furthermore, surge is easily extensible with more features and adaptable to further application.
Implementation
Data
We assembled a list of molecular formulae for benchmarking surge against MOLGEN 5.0 in Tables 1, 2. These formulae were taken from the natural products database COCONUT [22]. The size of these molecular formulae varies and is enough to challenge even the best constitutional isomer generators available (see Results section).
Algorithm and mathematical background
Surge is based on the nauty [23] package for computing automorphism groups of graphs as well as canonical labels. Like nauty, surge is written in a portable subset of C and runs on a considerable number of different systems.
Surge is an integration of three existing tools from the nauty suite [24]: (a) geng generates simple graphs based on certain boundary conditions, (b) vcolg colors vertices in the output of geng and (c) multig inserts multiedges in the output of the first two tools (Fig. 1).
An isomorphism between two graphs is a bijection between their vertex sets that maps edges onto edges. If the graphs have adornments, such as atom types for the vertices or bond multiplicities for the edges, then those adornments must be preserved by the mapping. If the two graphs are the same; i.e., the isomorphism is from a graph to itself, it is called an automorphism. The automorphisms form a group under the operation of function composition, called the automorphism group (Fig. 2).
The meanings of isomorphism and automorphism are different for each of the three stages in our algorithm. Referring to Fig. 1, at the first stage (which we call a simple graph) there are no vertex or edge adornments and all rotations and reflections, 10 in total, are automorphisms. When vertex adornments are added in the second stage, the atom type becomes significant so only the identity mapping and the reflection through the oxygen atom are automorphisms. In the final stage, edge adornments are added but in this example the automorphism group is not further reduced since the reflection through the oxygen atom preserves both atom type and bond multiplicity. Note how the automorphism groups in the second and third stages are subgroups of the automorphism groups in the previous stages.
First stage
Input to surge consists of a molecular formula such as C_{7}H_{12}O_{2}S. Based on the element counts, in this case C = 7, O = 2, S = 1, H = 12, the atom valencies are used to calculate the plausible range of the number of edges of a connected simple graph representing the topology of a molecule with this formula, with hydrogen atoms omitted. Then geng is called to generate all the connected simple graphs with those parameters, subject also to a maximum degree condition depending on the molecular formula [25]. Geng generates one graph from each isomorphism class and these are passed to the second stage as they are produced, without any need to store them [25]. In this example, there are 10 nonhydrogen atoms and the number of edges is in the range 9–11.
Second stage
Given a simple graph G from the first stage, the second stage assigns elements to vertices in all distinct ways. The element counts must be correct, and we must have valence \(\ge\) degree at each vertex. More onerously, we only want one member of each equivalence class of element assignment under the automorphism group of G (Fig. 3). We next explain how this is accomplished.
The vertices of G are arbitrarily numbered 1,2,…,n. An element assignment can be represented as a list showing the element assigned to each vertex in order of vertex number. For example, a valid list might be L = (C,C,C,S,O,C,C,C,O,C).
Automorphisms of G have an action on lists that permutes their entries. Namely, for list L and automorphism \(\gamma ,\) the list \(\gamma\)(L) assigns the same element to vertex \(\gamma\)(v) as L assigns to v, for each vertex v. Thus,
If L is a list of elements and \(\gamma\) is an automorphism, L and \(\gamma\)(L) give equivalent assignment of elements to the vertices of G. Our task in this stage is to choose exactly one assignment from each equivalence class. Given a fixed ordering of the elements, for example C < O < S, two lists can be compared lexicographically, for example
This enables us to define
the maximum list in the equivalence class of L. Note that canon(L) = canon(L’) if L and L’ are equivalent, so there is a unique maximum list L* in the equivalence class and we can recognize it by the condition canon(L*) = L*. To put it another way, if \(\gamma\)(L) > L for any automorphism \(\gamma\) then L \(\ne\) L*; otherwise L = L*.
Now we describe the conceptual method for the second stage. For given G:
This algorithm is efficient if the automorphism group Aut(G) is small, but that is not always the case. Therefore, we adopt a more complex approach. An automorphism of G is called minor if it merely swaps two leaves (vertices of degree 1) that have a common neighbour. The minor subgroup M \(\le\) Aut(G) is the subgroup generated by all the minor automorphisms.
A flower is a maximal set of leaves with the same neighbour. In the left graph of Fig. 4, the flowers are {1,2,3}, {6,10} and {9,11}. The minor subgroup M consists of all automorphisms that preserve the flowers, such as (1 2 3)(9 11). The order of M is \(3!\times 2!\times 2! = 24\). In addition to M, the automorphism group may contain automorphisms that do not preserve the flowers, such as (6 11)(7 8)(9 10). To capture such automorphisms, we colour the graph as in the right side of Fig. 4. Vertices not in flowers are coloured black. Within each flower, vertices are coloured red, blue, green, … in order of vertex number, using a fixed list of colours that does not include black. Now let N be the group of automorphisms that respect the vertex colours. In the example, N has only the identity and (6 9)(7 8)(10 11).
An arbitrary automorphism of G can be obtained by first applying an element of N to capture how the flowers are mapped to each other, and then applying an element of M to capture the movement of leaves within each flower. In both steps the choice is unique, so we have a factorization
(In the language of group theory, M is a normal subgroup and N is a complete set of coset representatives.) In the example, consider (1 2)(6 11)(7 8)(9 10). It swaps the flowers {6,10} and {9,11} so we choose the element of N which does that, namely \(\gamma\) = (6 9)(7 8)(10 11). Then we have to arrange the leaves within the flowers with an element of M, namely \(\delta\) =(1 2)(6 10)(9 11). This achieves \(\gamma \delta\) = (1 2)(6 11)(7 8)(9 10).
The main advantage of factoring Aut(G) = NM is the following.
Theorem
For any list L, L = canon(L) if and only if L = max { \(\delta\)(L)  \(\delta\) in M} and L = max { \(\gamma\)(L)  \(\gamma\) in N}.
Proof
The “only if” direction is obvious since M and N are subsets of Aut(G). Suppose in the other direction that L = max {\(\delta\)(L)  \(\delta\) in M} and L = max {\(\gamma\)(L)  \(\gamma\) in N}. From the factorization of Aut(G) we know that L* = \(\delta\)(\(\gamma\)(L)) for some \(\gamma\) in N and \(\delta\) in M. Note that in both L and L* the elements are in nonincreasing order within each flower, as they are maximized with respect to M. Also recall that the automorphisms in N preserve the order of vertex numbers within the flowers, by virtue of the fact that we coloured the vertices in order of vertex number when we computed N. This means that we can take \(\delta\) to be identity, and so L* = \(\gamma\)(L). This proves that L* = L, since L = max { \(\gamma\)(L)  \(\gamma\) in N}.
In order to implement the condition L = max { \(\gamma\)(L)  \(\gamma\) in M}, we don’t need to compute M explicitly. Instead, since M is generated by transpositions, it suffices that within each flower the elements are in decreasing order relative to vertex number. Using the ordering of elements that we have chosen, in the example we just need to enforce the inequalities element(1) ≥ element(2) ≥ element(3), element(6) ≥ element(10) and element(9) ≥ element(11). The program recursively assigns elements to vertices in order of vertex number and enforces these inequalities as they become active rather than at the end.
To implement the condition L = max { \(\gamma\)(L)  \(\gamma\) in N}, we compute N using nauty and test that \(\gamma\)(L) \(\le\) L for each \(\gamma\) in N. This is efficient in practice because N is very small most of the time.
We can also partly enforce N by means of inequalities: since vertex 6 is the least vertex in a nontrivial orbit {6, 9} of N, we can assume element(6) ≥ element(9). This is not necessary but it gives a small time improvement.
As an example, C_{7}H_{14}N_{2}O_{7} has 15,425,657,612 isomers. Using the factorisation Aut(G) = NM reduces the number of nontrivial groups processed by 58% and the maximum group size from 2592 to 72. The overall generation time is 18% less. In typical cases, the method provides about 10–40% reduction in cost.
Third stage
After the assignment of elements to vertices is complete, the program moves to the next stage of selecting a bond multiplicity for each edge. This is the same type of problem as in the second stage. Instead of a list of elements for each vertex, we have a list of multiplicities for each edge. Instead of Aut(G), we use the subgroup of Aut(G) that preserves the element assignment. Otherwise M and N are defined as before. In the implementation, we don’t use nauty to compute N but instead filter the N subgroup from the second stage, rejecting those automorphisms which don’t preserve elements and converting the others to their action on the edges.
The constraints we have at this time are that for each atom the total number of incident bonds counting multiplicity must be at most the valence of the atom, and that the total of (valence—incident bonds) over all atoms must equal the desired number of hydrogen atoms. Once these constraints are satisfied, there is exactly one way to add hydrogens (though the program does not add them explicitly).
As an example, geng makes 534,493 unlabelled simple graphs in 1.3 s for Lysopine C_{9}H_{18}N_{2}O_{4}. For these graphs, the second stage subgroup N is trivial 58% of the time and never larger than 72. Assignment of elements to vertices produces 3,012,069,151 vertexlabelled graphs in 90s.The N subgroup for the third stage is trivial 98% of the time and never larger than 24. Finally, the assignment of bond multiplicities produces 5,979,199,394 completed molecules in an additional 100 s.
As demonstrated by our examples, surge can generate molecular structures very quickly, allowing for the inspection of extremely large sets of isomers. The generation speed is several times faster than even the fastest output format (SMILES). On the other hand, any particular application will likely have stronger restrictions on the structure than just a molecular formula. For example, some substructures may make the molecule unstable or give it chemical properties undesirable in the application. Or, experimental investigation of an unknown compound may have determined some features of the structure, so that only molecules with those features are of interest.
For these reasons, surge provides a number of filters to limit the output. The 3stage generation method allows some of them to be implemented almost for free, and all of them are much more efficient than filtering the output through an external program. For example, restrictions on the number of short rings and the planarity of the molecule can be enforced at Stage 1. Surge also provides some "badlists" of forbidden substructures (many of them inspired by the corresponding feature of MOLGEN).
The opensource nature of surge allows for a more advanced feature. By writing small code snippets, the user can insert custom filters into any of the three stages, and also perform such tasks as adding extra elements and commandline options. Several worked examples are provided with the program.
Results
Surge is available under a liberal opensource License (Apache 2.0) on GitHub at https://structuregenerator.github.io as well as from https://users.cecs.anu.edu.au/~bdm/surge/.
The system can be built with the standard Unix Configure/Make scheme and the resulting standalone executable is then run from the command line. By default, surge generates all constitutional isomers of a given molecular formula. Surge can write output in either SDfile [26] or SMILES [27] format. SMILES output is produced very efficiently by constructing a template for each simple graph at the first stage, so that only atom types and bond multiplicity must be filled in before output.
We benchmarked surge with the set of molecular formulae given in Table 1. Since our motivation for developing structure generators is for the generation of large molecules, Table 1 consists of natural products, randomly selected from the natural products database COCONUT [22]. For the list of molecular formulae, surge outperformed MOLGEN by orders of magnitude (Fig. 5) and MOLGEN terminated at a builtin limit of 2^{31}–1 structures. Reported computation times were linearly extrapolated based on the MOLGEN timing for 2^{31}–1 structures and the actual number of isomers reported by surge. Note that surge generates between 7 and 22 million molecules per second for all of these examples.
Surge has a tiny memory footprint irrespective of the molecule size or the number of isomers. All of the examples in this paper run in at most 5 MB of RAM on Linux (Fig. 6).
For randomly selected 10 molecular formulae, 4 options of surge were tested and results are given in Table 2. These options are.

p0:1 At most one cycle of length 5

P The molecule is planar

B5 No atom has two double bonds and otherwise only hydrogen neighbours

B9 No atom lies on more than one cycle of length 3 or 4.
Limitations
Release 1.0 of surge does not perform a Hückel aromaticity test and therefore generates duplicate structures for Kekulé versions of aromatic rings that are graphtheoretically different. Benchmarking against MOLGEN 5.0 was therefore performed with the noaromaticity switch of MOLGEN.
Conclusion
We have presented surge, a structure generator for constitutional isomers based on the canonical generation path method. To the best of our knowledge, surge is the fastest chemical structure generator available. A number of badlist options are available to avoid the generation of potentially unlikely structures. Current limitations include the lack of an aromaticity detection. Surge is hosted as an opensource package on GitHub, inviting the scientific community to use and extend it. Surge offers a plugin mechanism for communitydriven extensions. Plugins can hook into the various stages of the surge generation process, thereby offering efficient means to prune the generation tree.
Availability of data and materials
Project name: surge
Project home page: https://structuregenerator.github.io
Operating system(s): Platform independent
Programming language: C
License: Apache 2.0
References
Elyashberg M, Argyropoulos D (2020) Computer assisted structure elucidation (CASE): current and future perspectives. Magn Reson Chem. https://doi.org/10.1002/mrc.5115
Miyao T, Kaneko H, Funatsu K (2016) Ring systembased chemical graph generation for de novo molecular design. J Comput Aided Mol Des 30:425–446
SaldívarGonzález FI, HuertaGarcía CS, MedinaFranco JL (2020) Chemoinformaticsbased enumeration of chemical libraries: a tutorial. J Cheminform 12:64
Blum LC, Reymond JL (2009) 970 Million druglike small molecules for virtual screening in the chemical universe database GDB13. J Am Chem Soc 131:8732–8733
Lindsay RK, Buchanan BG, Feigenbaum EA, Lederberg J (1993) DENDRAL: a case study of the first expert system for scientific hypothesis formation. Artif Intell 61:209–261
Gulyaeva KA, Artemieva IL (2020) The ontological approach in organic chemistry intelligent system development. Advances in Intelligent Systems and Computing. Springer, Singapore, pp 69–78
Badertscher M, Korytko A, Schulz KP, Madison M, Munk ME, Portmann P et al (2000) Assemble 2.0: a structure generator. Chemometrics Intellig Lab Syst. 51:73–79
Holt DF, Eick B, O’Brien EA (2005) Handbook of computational group theory. CRC Press, Boca Raton
Kreher DL, Stinson DR (2020) Combinatorial algorithms: generation, enumeration, and search. CRC Press, Boca Raton
Serov VV, Elyashberg ME, Gribov LA (1976) Mathematical synthesis and analysis of molecular structures. J Mol Struct 31:381–397
Molchanova MS, Shcherbukhin VV, Zefirov NS (1996) Computer generation of molecular structures by the SMOG program. J Chem Inf Comput Sci 36:888–899
Yirik MA, Steinbeck C (2021) Chemical graph generators. PLoS Comput Biol 17:e1008504
Faulon JL (1992) On using graphequivalent classes for the structure elucidation of large molecules. J Chem Inf Comput Sci 32:338–348
Faulon JL (1994) Stochastic generator of chemicalstructure. 1. Application to the structure elucidation of large molecules. J Chem Inf Comput Sci 34:1204–1218
Junker J (2011) Theoretical NMR correlations based structure discussion. J Cheminform 3:27
Nuzillard JM, Georges M (1991) Logic for structure determination. Tetrahedron 47:3655–3664
Gugisch R, Kerber A, Kohnert A, Laue R, Meringer M, Rücker C, et al. MOLGEN 5.0, a Molecular structure generator in advances in mathematical chemistry. Advances in mathematical chemistry; Basak, SC, Restrepo, G , Villaveces, JL, Eds.
Grund R, Kerber A, Laue R (1996) Construction of discrete structures, especially isomers. Discrete Appl Math 67:115–126
Grüner T, Laue R, Meringer M (1997) Algorithms for group actions: homomorphism principle and orderly generation applied to graphs. DIMACS Ser Discrete Math Theoret Comput Sci 28:113–122
Yirik MA, Sorokina M, Steinbeck C (2021) MAYGEN: an opensource chemical structure generator for constitutional isomers based on the orderly generation principle. J Cheminform. https://doi.org/10.1186/s13321021005299
Ruddigkeit L, van Deursen R, Blum LC, Reymond JL (2012) Enumeration of 166 billion organic small molecules in the chemical universe database GDB17. J Chem Inf Model 52:2864–2875
Sorokina M, Merseburger P, Rajan K, Yirik MA, Steinbeck C (2021) COCONUT online: collection of open natural products database. J Cheminform 13:2
McKay BD, Piperno A (2014) Practical graph isomorphism. II J Symb Comput 60:94–112
McKay B, Piperno A. nauty and Traces User’s Guide. 2019 Sep. https://pallini.di.uniroma1.it/Guide.html
McKay BD (1998) Isomorphfree exhaustive generation. J Algorithms 26:306–324
CTFILE FORMATS BIOVIA DATABASES 2016. 2016. https://help.accelrysonline.com/ulm/onelab/1.0/content/ulm_pdfs/direct/reference/ctfileformats2016.pdf
Weininger D (1988) SMILES, a chemical language and information system. 1. Introduction to methodology and encoding rules. J Chem Inf Comput Sci 28:31–36
Funding
Open Access funding enabled and organized by Projekt DEAL. MAY and CS acknowledge funding by the CarlZeissFoundation.
Author information
Authors and Affiliations
Contributions
BDM wrote the code and developed the underlying nauty package. BDM, CS and MAY conceived the project. BDM and CS guided the development. MAY contributed to the conceptual development and performed the evaluation and testing. All authors wrote the manuscript. All authors read and approved the final manuscript.
Corresponding authors
Ethics declarations
Competing interests
All authors declare no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
McKay, B.D., Yirik, M.A. & Steinbeck, C. Surge: a fast opensource chemical graph generator. J Cheminform 14, 24 (2022). https://doi.org/10.1186/s13321022006049
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s13321022006049
Keywords
 Structure generation
 Constitutional isomers
 Canonical generation path