Evolutionary Monte Carlo of QM Properties in Chemical Space: Electrolyte Design

Optimizing a target function over the space of organic molecules is an important problem appearing in many fields of applied science but also a very difficult one due to the vast number of possible molecular systems. We propose an evolutionary Monte Carlo algorithm for solving such problems which is capable of straightforwardly tuning both exploration and exploitation characteristics of an optimization procedure while retaining favorable properties of genetic algorithms. The method, dubbed MOSAiCS (Metropolis Optimization by Sampling Adaptively in Chemical Space), is tested on problems related to optimizing components of battery electrolytes, namely, minimizing solvation energy in water or maximizing dipole moment while enforcing a lower bound on the HOMO–LUMO gap; optimization was carried out over sets of molecular graphs inspired by QM9 and Electrolyte Genome Project (EGP) data sets. MOSAiCS reliably generated molecular candidates with good target quantity values, which were in most cases better than the ones found in QM9 or EGP. While the optimization results presented in this work sometimes required up to 106 QM calculations and were thus feasible only thanks to computationally efficient ab initio approximations of properties of interest, we discuss possible strategies for accelerating MOSAiCS using machine learning approaches.


Introduction
Increasing efficiency and longevity of energy storage systems is critical for improving economic sustainability of lowering greenhouse gas emissions. 1One aspect of this problem is searching chemical compound space for organic molecules optimal for a target application, such as lithium battery electrolyte component [2][3][4][5][6] or electroactive molecules for redox flow batteries. 7,8In this work we focused on the former, more specifically on finding electrochemically stable organic molecules that are good solvents for alkali salts.While such searches can be aided with high-throughput screening, [2][3][4] there has been a surge of ways to go beyond by increasing efficiency of compound property evaluations, e.g. with machine learning 9 or quantum alchemy, [10][11][12] and by sampling chemical space more efficiently.In the context of optimizing small organic molecules, most methods of the latter category can be classified as those based on Markov decision processes, [13][14][15][16][17][18] recurrent neural networks, 19,20 genetic algorithms, [21][22][23][24][25][26] and variational autoencoders. 27,28hile several variants of Markov chain Monte Carlo 29 sampling have also been applied to molecule optimization problems, 30,31 one intriguing variant, namely Evolutionary Monte Carlo, [32][33][34] has been overlooked so far.9][40] As illustrated in Figure 1, Evolutionary Monte Carlo involves running several Markov chain Monte Carlo simulations that focus on exploitation (i.e.refining already known molecules via incremental changes) or exploration (i.e.finding promising regions of chemical space), which interact by swapping configurations analogously to parallel tempering or by creating "child configurations" similarly to genetic algorithms in a way that observes detailed balance condition. 41s is the case for genetic algorithms, increasing the number of replicas yields more opportunities for creating "child configurations," thus accelerating exploration of chemical space.Unlike genetic algorithms though, Evolutionary Monte Carlo allows straightforward control of its exploration and exploitation aspects while guaranteeing to eventually find the global minimum due to the properties of Markov chain Monte Carlo.][23][24] While some recently proposed methods for molecular optimization operate in string representations, 18,19,25,52,53 we performed all procedures directly on chemical graphs to facilitate ensuring validity of generated molecules and maintaining detailed balance, as well as provide a more direct connection between the molecules considered and graph-based representations that have proven efficient in machine learning applications. 54,55Lastly, we implemented a simple Wang-Landau biasing potential 56 as a curiosity reward 57 increasing exploration aspect of the algorithm by "pushing" our Markov chain Monte Carlo simulations out of previously occupied graphs.The resulting method is named MOSAiCS (Metropolis Optimization by Sampling Adaptively in Chemical Space).][60] The rest of the paper is organized as follows.Section 2 presents the main ideas behind our approach in Subsecs.2.1-2.3,following up with description of the optimization problem on which we test it in Subsec.2.4 and details of our Monte Carlo simulations in Subsec.2.5.Section 3 discusses our experimental results, Section 4 concludes the paper with a results summary and outline of possible strategies to improve our approach.Some technical details of our method's implementation, experimental setup, and results were left for Supporting Information.
Figure 1: Scheme of an Evolutionary Monte Carlo workflow for molecular optimization that introduces several replicas with varying degree of focus on exploitation or exploration aspects of the calculation and evolves them by applying elementary mutations (blue) to single replicas and crossover moves (red) to pairs of replicas in a way satisfying detailed balance.

Chemical space definition
We aim to minimize a loss function F over a set of molecules, the latter represented by their chemical graphs.We define a chemical graph as an undirected graph whose nodes correspond to heavy atoms, along with, where present, hydrogen atoms covalently connected to them, and whose edges connect a pair of nodes if their heavy atoms share a covalent bond.For a chemical graph we also define a resonance structure as a set of valences of nodes' heavy atoms and orders of covalent bonds connecting these heavy atoms, both quantities taking integer values.The sum of covalent bond orders connecting a heavy atom to other atoms equals its valence, with the orders of bonds between a heavy atom and a hydrogen atom counted as one.Valence numbers are chosen to be chemically reasonable (e.g.IV for C, II, or IV, or VI for S) and we require their sum to be the minimum needed to build a set of covalent bond orders.We also forbid a covalent bond order to be larger than three.The reasons for not including valences and bond orders in the definition of a chemi-cal graph, but rather enumerating their possible values separately, are illustrated in Figure 2 demonstrating examples of molecules for which several resonance structures differing in bond orders or bond orders and heavy atom valences can be defined.We emphasize that while this definition of bond orders and valences is loosely based on valence structure theory, it was designed not to reflect actual electronic structure of a molecule, but to allow convenient definitions of changes of chemical graphs that are illustrated in Figures 3 and 4, as will be discussed in detail in Subsec.2.3.
Our definition of chemical graph a priori prevents us from differentiating between conformers or stereoisomers and we will assume our optimization problem to be unaffected by this, e.g. if for a given chemical graph we are interested only in the most stable stereoisomer and we optimize a Boltzmann average.We also did not implement support for molecules where valid Lewis structures can only be generated by assigning charges to atoms, e.g.compounds with nitro groups, hence they were ignored during all calculations done in this work.

Monte Carlo sampling
We perform optimization by running a Markov chain Monte Carlo simulation (referred to as just "simulation" from now on) of unnormalized probability density similar to the one used for parallel tempering where X is a set of N repl chemical graphs (also referred to as replicas) . ., N repl ) are temperature parameters, V bias is biasing potential, C ∞ is the "infinitely convex function" defined to be such that for arbitrary sets of numbers x (j) and y (j) (j = 1, . . ., N opt ) min j=1,...,Nopt and N opt is the number of replicas that, as will become clear later, effectively undergo greedy stochastic minimization and are referred to as greedy replicas, with the other replicas, referred to as exploration replicas, providing a less restricted exploration of chemical space and preventing greedy replicas from getting stuck in a local minimum of F .The history-dependent biasing potential bias is defined as 56 where ρ (j) (X) is the number of times X has been visited during the simulation by replica with index j (details on how it was evaluated are left for Subsec.2.5), α b is the user-defined bias proportionality coefficient.Setting a nonzero α b makes sampling X non-Markovian; as a result our certainty that in this regime a global minimum of F w.r.t.X is eventually found is based not on properties of Markov chain Monte Carlo, but on heuristic expectation that the biasing potential would make probability distribution of each exploration replica approach uniformity, leading to at least one replica coming across the global minimum over a finite number of simulation steps.

Monte Carlo moves
A simulation consists of taking a sequence of moves in a way outlined in Algorithm 1.If the current set of replicas is in configuration X 1 , a move involves randomly generating parameters w of a change and deciding to replace X 1 with the change's outcome (or trial configuration) X 2 with an acceptance probability similar to the standard Metropolis-Hastings expression 41 where P prop (X 1 , w) is the probability that w is proposed given that X 1 is the initial configuration and w −1 are parameters of a random change yielding X 1 when applied to X 2 and corresponding to a unique w.The latter property ensures that detailed balance still holds in situations when several w yield the same trial configuration X 2 .Trial configurations decreasing the minimal value of F among greedy replicas compared to initial configurations is accepted automatically due to our definition of C ∞ (2).
Require: Initial configuration: X 1 ; loop Randomly choose change parameters w; Use w on X 1 to generate X 2 ; Randomly sample r from uniform distribution in [0, 1]; r acc.← P acc.(X 1 , w, X 2 ) (see Eq. ( 4); if r < r acc.then X 1 ← X 2 ; end if end loop We use three types of moves to propose the trial configurations X 2 ; we will only discuss the general idea behind them here with implementation details left for Supporting Information.The first type, referred to as elementary moves, applies an elementary mutation outlined in Figure 3 to a single replica; such moves correspond to incremental exploration of chemical space.To accelerate greedy optimization of molecules, we additionally introduced the "no reconsiderations condition": if change parameters w corresponding to an elementary move have been rejected for a greedy replica they are not considered again.The second type of moves is tempering swap moves that are analogous to the swap moves in conventional parallel tempering techniques and involve randomly choosing replicas with indices i and j in such a way that at least one of them is an exploration replica, considering a swap of the corresponding chemical graphs, and accepting it with acceptance probability (4).These moves allow greedy replicas stuck in a local minimum of F to get to chemical graphs with lower values of F if the latter are discovered by an exploration replica.
The third type of moves are crossover moves inspired by the procedure developed in Ref.21, which are introduced to allow drastic changes of chemical graphs occupied by replicas.The general idea is illustrated in Figure 4: a pair of nodes is randomly chosen in two chemical graphs and the neighborhoods of these two nodes are exchanged to create two new chemical graphs.Thus defined crossover moves are more restrictive than the ones of Ref. 21 as they do not allow exchanging fragments of arbitrary shape and connectivity.These restrictions, however, make it straightforward to ensure that the resulting chemical graphs satisfy constraints on the number of nodes, are connected, and correspond to a change for which the P prop (X 1 , w)/P prop (X 2 , w −1 ) ratio in P acc (4) can be easily calculated.
Setting F to be infinitely large for chemical graphs violating certain user-defined constraints is a general way to enforce the latter on the optimization result.However, it is in general preferable to maintain a given constraint as early as during the proposition of trial configuration X 2 to increase average acceptance probability and the resulting speed of chemical space exploration.We implemented the corresponding algorithms for maintaining constraints on the number of heavy atoms in a molecule and the kinds of atoms that can share a covalent bond since they are simple to maintain, yet quite important for our applications.Lastly, the question of the moves' sufficiency to access the chemical space and sets of molecules considered in Section 3 in their entirety is discussed in Supporting Information.

Minimization problems
A good battery electrolyte is a good solvent for lithium salts and is electrochemically stable.We approximated the former property with polarity; maximizing a molecule's polarity was in turn approximated by either maximizing the dipole moment D or minimizing the free energy of solvation in water ∆G solv .We approximated the electrochemical stability requirement with a lower bound on the HOMO-LUMO ∆ϵ, with which we approximated the width of the compound's electrochemical stability window. 2 While the latter relation is not actually practical for battery design, 61 we still opted for a ∆ϵ-based electrochemical stability criterion to M1: Addition/removal of a node (red) connected to a single neighbor (blue).
M2: Changing order of a covalent bond (green) and number of hydrogen atoms in the nodes it connects (blue).
M3: Replacing a heavy atom (red) in a node without changing its connections to other nodes.
M4: Change valence of a heavy atom (green) by adding or removing hydrogen atoms connected to it.
M5: Change valence of a heavy atom (green) by adding or removing nodes (red) connected to it.
M6: Change valence of a heavy atom (green) by changing order of a covalent bond (green) with another node whose hydrogen number is changed (blue).II -selecting the "blue" neighborhoods of the two blue nodes, coloring the rest "red," with the bonds connecting blue and red fragments colored "green"; III -exchanging the blue fragments between molecules and reconnecting them with red fragments by exchanging connected nodes between pairs of green bonds.
connect our work with other compound optimization problems where ∆ϵ can be used. 62,63or both D and ∆G solv optimization we constrained the molecules' ∆ϵ to be larger than either benzene (strong ∆ϵ constraint) or octa-1,3,5,7-tetraene (weak ∆ϵ constraint), resulting in four minimization problems of differing difficulty.While in this work we focused on testing performance of MOSAiCS against these single objective optimization problems, our approach can also be used to optimize several properties at once via a suitable multiobjective loss function. 64e aimed to estimate ∆G solv , D, and ∆ϵ as computationally cheaply as possible while being qualitatively correct over a wide range of chemical compounds; the resulting protocol is explained in detail in Supporting Information.Here we just mention that for a given chemical graph we used the MMFF94 forcefield 65 to generate molecular conformers, for which we performed GFN2-xTB 66 calculations with analytical linearized Poisson-Boltzmann model 67 simulating presence of water.The root mean square error (RMSE) that is presented for cal-culated quantities corresponds to statistical error from randomness of conformer generation.We used two sets of parameters for our protocol: "converged" that produced reasonable RMSEs for a wide variety of compounds, but was relatively computationally expensive, and "cheap" that was used during our simulations.From now on, ∆G solv , D, and ∆ϵ will denote estimates of these quantities obtained with the "converged" protocol, while estimates obtained with the "cheap" protocol will be marked with addition of "cheap" superscript.
Each of the four minimization problems was solved in two sets of molecules based on QM9 68 and the Electrolyte Genome Project 5 (EGP) datasets.The QM9 dataset consists of 134k molecules containing up to 9 heavy atoms (C, O, N, and F).We defined the "QM9*" set to consist of molecules (not necessarily in QM9) that also contain up to 9 heavy atoms of the same elements as QM9, but are additionally constrained by not allowing bonds between N, O, and F atoms, as well as O-H and H-F bonds, since these covalent bonds are typically associated with increased chemical reactivity.The EGP dataset was generated with the Materials Project 69 workflows in an effort to facilitate discovery of novel battery electrolyte molecules; the version currently hosted on the Materials Project website contains 24.5k species in total; neutral species for which MMFF94 coordinates could be generated included 19.7k individual chemical graphs containing up to 92 heavy atoms.These characteristics of the EGP dataset were the basis for defining the "EGP*" set, whose molecules (not necessarily in EGP) contain up to 15 heavy atoms (B, C, N, O, F, Si, P, S, Cl, and Br, which are elements present in organic molecules of EGP) and, for the sake of chemical stability, do not contain covalent bonds between N, O, F, Cl, and Br, between H and B, O, F, Si, P, S, Cl, or Br, as well as S-S and P-P bonds.
We chose 15 as the maximum number of heavy atoms allowed in EGP* molecules because this size restriction is obeyed by 87.0% and 97.0% of EGP's chemical graphs satisfying weak and strong ∆ϵ constraints.When choosing which elements can not share a covalent bond in QM9* and EGP* molecules we mainly aimed for excluding weak bonds, although we also forbade some relatively strong bonds whose presence can signify molecular reactivity.Since we only consider molecules whose valid Lewis structures can be generated without assigning charges to atoms, H-F and double O-O bonds can only be encountered in hydrogen fluoride and oxygen, which we excluded from consideration due to their corrosive properties.Creating N-N bonds inside an organic compound risks making it prone to releasing nitrogen on excitation, adding functional groups containing double N-O bonds to a molecule risks making the latter prone to self-oxidation, and hydroxyl groups engage relatively easily in reactions involving oxidation or nucleophilic attacks. 70We note that in practice, managing this kind of reactive behavior would require additional use of more sophisticated compound stability measures.
While both QM9* and EGP* are well defined and finite sets of chemical graphs, their huge size makes evaluating any of their properties exactly, i.e. through exact enumeration of all their chemical graphs, unfeasible.However, we do summarize properties of intersections of QM9* and QM9, as well as EGP* and EGP, in Supporting Information.

Simulation details
During a simulation we used ∆ϵ cheap to estimate whether a molecule satisfies the constraint on ∆ϵ; dimensionless loss functions corresponding to D and ∆G solv were defined as where STD dataset refers to standard deviation of a quantity over molecules at the intersection of chemical graph set of interest and the reference dataset (QM9 for QM9* and EGP for EGP*) which satisfy the ∆ϵ constraint of interest.We chose 1000 "pre-final" molecules exhibiting the smallest value of loss function out of the molecules visited during the simulation and evaluated converged estimates of the quantities of interest for them; the molecule with the best D or ∆G solv value among pre-final molecules satisfying the ∆ϵ constraint is the one considered the candidate molecule proposed by the simulation.We used N repl = 36 with N opt = 4 (cf.definitions in Subsection 2.2); virtual temperature parameters β (i) appearing in P (1) were defined in such a way that the smallest and largest β (i) were 1 and 8, and the other β (i) formed a geometric progression between the two extrema values, the latter being a simple recipe taken from applications of parallel tempering to configuration space sampling. 71,72A simulation consisted of 50000 "global" steps, out of which 60% were "simple" steps applying an elementary move to each replica, 20% were "tempering" steps making tempering swap moves on 128 randomly chosen pairs of replicas, and another 20% were "crossover" steps making crossover moves on 32 randomly chosen pairs of replicas.ρ (j) (X) appearing in V (i) bias (3)  was counted as the number of times replica j was found in X after a global step had been completed.For elementary moves we additionally set that: the nodes added or removed during M1 mutation could be connected to the molecule with bonds of order from 1 to 3; bonds changed with M2 and M6 mutations could have their order increased or decreased by 1 and 2 respectively; nodes added or removed with M5 mutation could be connected to the molecule with bonds of order 1 or 2.
We set α b to 0.0, 0.2, or 0.4; for each of the resulting 12 combinations of α b and optimization problem we ran 8 simulations with different random number generator seeds.For all simulations all replicas initially occupied the chemical graph of methane.While it would be natural to assign each replica a randomly chosen molecule from the intersection of QM9 and QM9* or EGP and EGP*, we went with the intentional handicap of using methane as the starting molecule to demonstrate that MO-SAiCS is capable of constructing all the candidate molecules presented in this Section from scratch.The effect of choice of initial condi-tions on the final result is briefly addressed in Supporting Information.

Results and discussion
In this Section we describe the main results of our numerical experiments.The more technical aspects, such as full information on generated candidates and influence of biasing potential on search efficiency, are left for Supporting Information.
While we ran in total 96 simulations in QM9*, or 24 simulations with different random generator seed and α b values for each optimization problem, they agreed remarkably often on candidates proposed, only yielding 10 candidates in total.Table 1 summarizes the best and worst values of optimized quantities of candidates proposed by MOSAiCS along with the corresponding relative improvement, which we define as absolute difference between a candidate's optimized quantity value and the corresponding value for the best molecule for the optimization problem taken from the reference dataset (cf.Table S2 in Supporting Information), divided by the corresponding STD dataset .For optimization of ∆G solv with weak ∆ϵ constraint all trajectories proposed the minimum of ∆G solv already present in QM9, while for all other optimization problems all trajectories proposed candidates that improved significantly on molecules in QM9.Best candidates proposed for a given optimization problem are shown in Figure 5; note that to facilitate discussion of candidates' properties in Supporting Information, each candidate is referred by a capital C with a unique index subscript and a superscript denoting the reference dataset.We see how MOSAiCS successfully constructed complex conjugated bond structures facilitating charge transfer which, in turn, led to smaller, i.e. more negative, ∆G solv or larger D. Figure 6 displays optimization progress with number of global Monte Carlo steps for different α b values for minimizing ∆G solv with weak ∆ϵ constraint.We observe convergence of the optimized property with a rate not significantly affected by changing α b ; the same is true with varying degree for other optimization problems as discussed in Supporting Information.To visualize how simulations explored chemical space for different optimization problems and values of α b , for each such combination we took a simulation that had produced the best candidate and plotted the density of molecules it encountered with respect to the optimized quantity and ∆ϵ. Figure 7 presents such plots for minimizing ∆G solv with weak ∆ϵ constraint, with plots for other optimization problems presented and discussed in Supporting Information.We see that increasing α b tended to increase diversity of molecules encountered during the simulations, but this was mainly done by considering more molecules in regions of chemical space with larger values of ∆G solv .Optimization in EGP* was harder than in QM9* due to larger size of the former set of molecules, resulting in our protocol generating underconverged simulations that rarely agreed on candidates, producing 83 candidates in total.However, as summarized in Table 2, we still observed significant improvement of optimized quantities compared to EGP, though the improvements' impressive magnitudes are largely due to EGP containing a much less representative portion of EGP* compared to the case  8; unlike QM9* candidates, no chemical intuition is seen in how they were constructed beyond adding as many polar covalent bonds as possible, which may be due to underconvergence of our EGP* simulations.The underconvergence can also be observed on our optimization progress plots, the plot for minimizing ∆G solv with weak ∆ϵ constraint presented in Figure 9 and the rest found in Supporting Information.Figure 9 also demonstrates how adding α b can accelerate optimization as a function of global Monte Carlo steps, though we need to note that simulations with larger α b on average process more chemical graphs per global Monte Carlo steps as discussed in Supporting Information.Densities of molecules encountered during simulations minimizing ∆G solv with weak ∆ϵ constraint that produced the best candidates are presented in Figure 10; unlike the case of QM9*, increasing α b helped simulations explore parts of chemical space with smaller values of ∆G solv .Analogous plots for other optimization problems in EGP* are presented in Supporting Information.

Conclusions and outlook
We have proposed an effective algorithm for optimization in chemical space, dubbed MO-SAiCS, and successfully applied it to several test optimization problems connected to lithium battery electrolyte design.In the cur-    rent implementation, it is only feasible to optimize estimates of quantities of interest that can be evaluated with little computational effort due to the large number of evaluations of loss function made during a simulation (see Supporting Information); given successes in using active learning for optimization problems in both configuration [73][74][75][76] and chemical 58,77,78 space, our first priority is to combine MOSAiCS with a similar protocol to decrease the number of loss function evaluations done during the simulations.Successful use of Markov Decision Process formalism to accelerate genetic algorithms in chemical space 59 suggests MOSAiCS might similarly be improved with a smarter policy for choosing elementary mutations and crossover moves.On a more general note, any method for generating chemical graphs that can also provide corresponding proposition probability P prop needed for P acc (4) can be integrated into MOSAiCS framework directly.
While we aimed to propose an approach that would be agnostic to how much is known about the chemical graph set of interest, we still relied on QM9 and EGP to get reasonably rescaled loss functions F solv. (5) and F dipole (6) that were then used during simulations in QM9* and EGP*.This dependence on previously published data should become avoidable by implementing more sophisticated schemes 79 for adjusting temperature parameters β (i) based on trajectory history.Also, while we used heavy atoms with connected hydrogens as nodes of chemical graphs to maximize chemical diversity of generated compounds, it is possible to expand the algorithm to using larger compound fragments as nodes instead.If applicable to the optimization problem at hand, this modification should simplify the search by both decreasing effective dimensionality of the graphs considered 26,31 and improving scalability of machine learning models for the molecules of interest.

TOC Graphic 1 Details of elementary mutation and crossover move implementations
While in Subsec.2.3 of the main text it was convenient to speak about invertible elementary mutations, the corresponding procedures were implemented in terms of elementary changes listed in the left column of Table S1, which differ from elementary mutations in treating procedures for creating and destroying nodes separately rather than part of unified mutations changing the number of nodes in a chemical graph.To apply a random mutation to a chemical graph we enumerate all sets of possible elementary change parameters, which is also when we observe constraints on allowed types of covalent bonds and number of nodes.
We then randomly choose change parameters among those enumerated, by first choosing the elementary change, and then all other parameters in the order listed in the right column of Table S1.Lastly, we check that straightforward application of the change yields a molecule without nodes with invalid valences and if that is not the case reject the move automatically.
The last step's necessity is illustrated with an example in Figure S1.
In Table S1, whenever we mention choosing a node or a pair of nodes we mean choosing them among a list of all nodes or pairs of nodes such that none of them are equivalent to any other in the chemical graph.This prevents MOSAiCS from considering several mutations that yield the same chemical graph X 2 given an initial chemical graph X 1 , making it straightforward to define the unique inverse change parameters that yield X 1 when applied to X 2 .Enumerating all possible changes for X 2 analogously to X 1 allows straightforward calculation of the proposition probability ratio in acceptance probability expression [Eq.(4)   in the main text].Choosing resonance structure for E2 and E6 is necessary in situations when decreasing bond order, depending on the resonance structure, can yield two distinct molecules corresponding to cases when the changed covalent bond is conserved or broken; such situations are illustrated in Figure S2.
Choosing the parameters of a crossover move starts with assigning the "blue fragments" S-2 Table S1: Elementary changes used in this work along with the corresponding order in which random change parameters were chosen.

Elementary change
Choice order add node (E1a)  (as illustrated in Figure 4) which are defined as the set of nodes separated by a given number of covalent bonds from an "origin" node.Origin nodes in the two molecules are chosen randomly; then we enumerate all possible blue fragments such that: • a pair of resonance structures chosen for both molecules allows a one-on-one mapping between green (see Figure 4) covalent bonds of the same order; • exchanging the blue fragments yields chemical graphs satisfying constraints on the number of nodes in a chemical graph; • at least one blue fragment and at least one red (see Figure 4) fragment contain more than one node, ensuring that the crossover move is not a simple combination of two E3 elementary changes; • the number of green bonds is not larger than a certain threshold (set 3 in this work); this is done to limit the number of chemical graph pairs enumerated during the next step.
With the blue fragments chosen we enumerate all pairs of chemical graphs obtained by exchanging them, i.e. taking pairs of green bonds from molecules and exchanging nodes between them, thus connecting each molecule's blue fragment to the red fragment of the other molecule.For each thus generated pair we additionally check that the valences in the exchanged blue and red fragments are valid in the new molecules (to prevent situations analogous to the one illustrated in Figure S1), that none of the newly created green bonds connect the same pair of nodes (to prevent irreversible changes illustrated in Figure S3), and that none of the newly created bonds violate constraints on which types of heavy atoms can be covalently connected.
Of the resulting pairs of chemical graphs one pair, consisting of chemical graphs X′ and X′′ , is randomly chosen.The chemical graphs are assigned to replicas i ′ and i ′′ according to the value of sampled probability P [Eq.( 1) in the main text] corresponding to the ordering.exhibited low RMSE's over QM9 and EGP datasets, but sometimes required relatively large computational time to be calculated, hence the cheap set of parameters were used during our simulations.
Properties of ∆G solv , D, and ∆ϵ over molecules in intersections of QM9 with QM9* and EGP with EGP* satisfying ∆ϵ constraints are summarized in Table S2.Note that RMSE's of properties of interest are reasonable over all such molecules.The table also presents reference values used to calculate relative improvements presented in the main text for candidate molecules.

Accessibility of chemical space
To guarantee (at least for α b = 0 when the sampling is Markovian) that MOSAiCS eventually finds the minimum of loss function over the set of chemical graphs of interest the procedures used to propose trial replica configurations should have connectivity, i.e. it should be possible to morph each pair of members of the set from one into another using these procedures alone and with only other molecules in this set as intermediates.Mutations M1-M4 are enough to guarantee connectivity in a large variety of relevant sets of chemical graphs, for example those obeying an upper bound on the number of heavy atoms and a constraint on which types of non-carbon heavy atoms can share a covalent bond.Connectivity of this kind of sets is easily seen by "constructing" and "deconstructing" each set member to and from methane.The same argument holds for QM9* but, unfortunately, not for EGP* due to P and S atoms being forbidden to share a covalent bond with a hydrogen atom, making it impossible to use mutation M4.Introducing additional mutations M5 and M6 alleviates the problem only partially, with an example of an EGP* molecule that cannot be transformed to and from methane using elementary mutations presented in Figure S4.This molecule can be considered by the simulations used in this work thanks to crossover moves, but we cannot guarantee this is the case for all such examples.However, we note that all such molecules are geometrically strained, which makes them not interesting for our applications due to S-8 Table S2: Properties of sets of molecules at the intersection of QM9 and QM9* or EGP and EGP* that satisfy weak or strong constraint on HOMO-LUMO gap ∆ϵ, namely their number (num.mol.), standard deviations (STD) of optimized quantities (free energy of solvation ∆G solv and dipole D) over these molecules, the maximum root mean square errors (RMSEs) observed ∆G solv , D, and ∆ϵ, minimum of ∆G solv and maximum of D, along with SMILES of the molecules for which these extrema values were observed.3 Additional data obtained during optimization in QM9* and EGP* Full information about QM9* candidates, including the corresponding values of optimized quantities and how often they were proposed by our simulations, is presented in Table S3; we see that for all optimization problems but for maximizing D with strong ∆ϵ constraint the agreement between simulations in terms of proposed candidates is close to unanimous.For QM9* simulations mean and standard deviation values of these quantities are summarized in Table S4.Changing the bias proportionality coefficient α b did not significantly affect the number of global steps required on average to find a candidate molecule, but did lead to S-10 a more thorough (but also costly) exploration of chemical space, as demonstrated by significant increase of N graph req and N graph tot .Another interesting observation is that for a given ∆ϵ constraint N step req values are significantly higher for optimization of D than optimization of ∆G solv .This is likely due to ∆G solv being easier to optimize since free energy of solvation can be accurately decomposed as the sum of scalar contributions from different functional groups in the molecule, S11 while for the dipole moment such contributions would take vector form whose sum would depend on the groups' placement in Cartesian space.Also note that for all optimization problems but for maximizing D with a strong ∆ϵ constraint the N step req was typically significantly smaller than the simulation length, which is an indicator of the simulations being converged in these cases.Relative improvement progress plots for optimization problems not shown in the main text are presented in Figure S6; while they look largely analogous to Figure 6 (1.523 ± 0.000) (1.036 ± 0.007) (1.000 ± 0.000) • 10 1 2 0 0     S5; we observe from the behavior of average optimized property values of candidates D best and ∆G best solv that adding biasing potential at worst did not significantly affect optimization results and at best provided significant improvements.Lastly, the difference between cheap and converged estimates of D and ∆G solv is very small for rigid molecules that only have one local minimum in configuration space (e.g. benzene), but it can become significant for larger non-rigid molecules with many valid conformers.To compare how much this difference affected candidate molecules obtained in this work we introduced "numerical error measure" ∆ conv cheap defined as the ratio of mean absolute error of the cheap estimate relative to the converged one divided by standard deviation of the converged estimate over the pre-final (as defined in Subsec.2.5 of the main text) molecules.
Its values for QM9* and EGP* simulations are presented in Table S6.For QM9*, we note that maximizing D under strong ∆ϵ constraint was apparently much more affected by using D cheap in the minimization function, which may be the reason simulations for that optimiza-S-16

S-18
tion problem agreed on candidate molecules much less frequently than for the other three.
We also observe that for a given optimization problem simulations in EGP* were much more affected by cheap/converged difference than the ones in QM9*, likely due to EGP* consisting of larger molecules with more conformers to be accounted for.This added to EGP* being the more challenging set of molecules to search through compared to QM9*.

Effect of initial replica positions on optimization results
It is natural to assume that increasing diversity of molecules initially occupied by replicas should make it easier for MOSAiCS to propose better candidate molecules.We tested this hypothesis on minimizing ∆G solv with strong ∆ϵ constraint over QM9*, since it is a problem for which converged results proved easy to obtain and for which improvement over QM9 was observed.For that optimization problem we launched 8 simulations with α b = 0.0 which only differed from the others presented in this work by the choice of initial molecules: each was chosen with uniform probability from molecules in the intersection of QM9 and QM9* that S-19 satisfied strong constraint on ∆ϵ cheap , then initial molecules were assigned to replicas in such a way that molecules with larger ∆G cheap solv were occupied by replicas with smaller β (i) .The resulting mean optimization progress is plotted in Figure S14, with the corresponding values observed when all replicas were initialized in methane added for comparison.Surprisingly, while simulations with random starting molecules had significantly larger starting values of relative improvement, they failed to go beyond those values significantly earlier than simulations initialized with methane, and they all proposed previously mentioned C QM9 2 as the candidate.While there might be better ways to choose initial molecules that would make search for candidates tangibly faster, finding them was beyond the scope of this work.

Figure 2 :
Figure 2: Examples of molecules that can be written in terms of several resonance structures that differ in (a) covalent bond orders and (b) both covalent bond orders and heavy atom valences.

Figure 3 :Figure 4 :
Figure 3: Definitions and examples of elementary mutations.A node is colored green if its heavy atom changes valence, otherwise it is colored red if it is destroyed or created, or colored blue if it changes the number of hydrogen atoms connected to the heavy atom.A bond is colored red if it is destroyed or created or green if it, in general, only changes bond order, though the latter may involve changing the order to and from 0. The valences in M5 and M6 always change to adjacent values (e.g. for S valence can change from II to IV, but not from II to VI); nodes created and destroyed in M1, M3, and M5 contain heavy atoms in their smallest valence state (e.g.valence II for S).

6 Figure 5 :
Figure 5: QM9* candidates that exhibited smallest free energy of solvation ∆G solv or largest dipole D under weak or strong constraint on the HOMO-LUMO gap ∆ϵ.See Supporting Information for more information about them.

Figure 7 :
Figure 7: Densities of molecules encountered by example simulations minimizing free energy of solvation ∆G solv with weak HOMO-LUMO gap ∆ϵ constraint in QM9* for different values of bias proportionality coefficient α b .

Figure 8 :
Figure 8: EGP* candidates that exhibited smallest free energy of solvation ∆G solv or largest dipole D under weak or strong constraint on the HOMO-LUMO gap ∆ϵ.See the Supporting Information for more information about them.

Figure 9 :
Figure 9: Relative improvement observed at N step global Monte Carlo steps for EGP* simulations optimizing ∆G solv with weak ∆ϵ constraint; results are organized analogously to Figure 6.

Figure 10 :
Figure 10: Densities of molecules encountered by simulations minimizing free energy of solvation ∆G solv with weak HOMO-LUMO gap ∆ϵ constraint in EGP* that produced the candidate with smallest ∆G solv for a given value of bias proportionality coefficient α b .
Details of MOSAiCS implementation and quantity of interest evaluation, discussion of accessibility of QM9* and EGP* by our simulations, additional experimental results.Acknowledgement This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 957189 (BIG-MAP) and No 957213 (BATTERY 2030+).O.A.v.L. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 772834).O.A.v.L. has received support as the Ed Clark Chair of Advanced Materials and as a Canada CIFAR AI Chair.O.A.v.L. acknowledges that this research is part of the University of Toronto's Acceleration Consortium, which receives funding from the Canada First Research Excellence Fund (CFREF).Ob-taining the presented computational results has been facilitated using the queueing system implemented at http://leruli.com.The project has been supported by the Swedish Research Council (Vetenskapsrådet), and the Swedish National Strategic e-Science program eSSENCE as well as by computing resources from the Swedish National Infrastructure for Computing (SNIC/NAISS).

Figure S1 :Figure S2 :
Figure S1: Example of an elementary change that is irreversible due to our definition of valid valence of an atom.Applying elementary change E2 to chemical graph A yields chemical graph B, which is actually an incorrectly written chemical graph C. Applying the inverse of A→B to C yields D, which differs from A.

Figure S4 :
Figure S4: An EGP* molecule that cannot be transformed to or from methane through EGP* molecules using elementary mutations alone.

Figure 7 , C QM9 8 , and C QM9 9 demonstrate
Figure S5 presents QM9* candidates that were not plotted in Figure 5; while we once again see MOSAiCS' ability to generate unconventional molecular structures, candidates C QM9 7 ,

10 Figure S5 :
Figure S5: QM9* candidates that did not exhibit best target property values.

Figure S6 :
Figure S6: Relative improvement (as defined in Sec. 3 of the main text and estimated from ∆G cheap solv or D cheap ) observed at N step global Monte Carlo steps for QM9* simulations optimizing (a) D with weak ∆ϵ constraint, (b) ∆G solv with strong ∆ϵ constraint, and (c) D with strong ∆ϵ constraint.For each value of the biasing potential α b we plot the mean value over different random generator seeds, the error bars corresponding to standard deviation."init.val." is the relative improvement of the simulations' starting molecule, i.e. methane.

Figure S7 :
Figure S7: Densities of molecules encountered by example simulations minimizing free energy of solvation ∆G solv with strong HOMO-LUMO gap ∆ϵ constraint in QM9* that produced the candidate with smallest ∆G solv for a given value of bias proportionality coefficient α b .

Figure S8 : 14 Figure S9 :
Figure S8: Densities of molecules encountered by example simulations maximizing dipole D with weak HOMO-LUMO gap ∆ϵ constraint in QM9* that produced the candidate with largest D for a given value of bias proportionality coefficient α b .

Figure S10 :
Figure S10: Relative improvement observed at N step global Monte Carlo steps for EGP* simulations optimizing (a) D with weak ∆ϵ constraint, (b) ∆G solv with strong ∆ϵ constraint, and (c) D with strong ∆ϵ constraint.Results are organized analogously to Figure S6.

Figure S11 :
Figure S11: Densities of molecules encountered by simulations minimizing free energy of solvation ∆G solv with strong HOMO-LUMO gap ∆ϵ constraint in EGP* that produced the candidate with smallest ∆G solv for a given value of bias proportionality coefficient α b .

Figure S13 :
Figure S13: Densities of molecules encountered by simulations maximizing dipole D with strong HOMO-LUMO gap ∆ϵ constraint in EGP* that produced the candidate with largest D for a given value of bias proportionality coefficient α b .

Figure S14 :
Figure S14: Relative improvement observed at N step global Monte Carlo steps for QM9* simulations optimizing ∆G solv with strong ∆ϵ constraint if the molecules initially occupied by replicas are chosen randomly from QM9 dataset ("rand.init.mol.") with bias proportionality coefficient α b set to 0.0.Plotted for comparison are corresponding values when all replicas are initialized at methane ("CH 4 init.mol.", also appearing in Figure S6a) as well as the relative improvement value for methane ("CH 4 val.").

Table 1 :
Best and worst QM9* candidates proposed during minimization of free energy of solvation ∆G solv or maximization of dipole D with weak or strong constraint on the HOMO-LUMO gap ∆ϵ, along with their optimized quantity values and the relative improvement compared to QM9 dataset as defined in Sec. 3. ∆G solv and D values are in kJ/mol and debye.The full list of candidate molecules can be found in Supporting Information.step global Monte Carlo steps for QM9* simulations optimizing ∆G solv with weak ∆ϵ constraint.For each bias proportionality coefficient α b we plot the mean over different random generator seeds, the error bars corresponding to standard deviation."init.val." is the relative improvement of the simulations' starting molecule, i.e. methane.
Figure 6: Relative improvement (as defined in Subsec.3 and estimated from ∆G cheap solv ) observed at N

Table 2 :
Best and worst EGP* candidates, data notation is analogous to Table1.Full list of candidate molecules can be found in Supporting Information.
The total number of molecules (tot.num.mol.) is the number of molecules whose chemical graphs, if extractable from QM9 or EGP coordinates with xyz2mol code, satisfied constraints of QM9* or EGP* sets.

Table S6 :
Measure of error caused by using cheap version of optimized quantity during the simulation ∆ conv cheap (see Sec.

Table S7 :
EGP* candidates proposed during minimization of free energy of solvation ∆G solv with weak constraint on the HOMO-LUMO gap ∆ϵ, the corresponding values of ∆G solv , and the number of times a molecule was proposed as a candidate by simulations run with different values of the bias proportionality coefficient α b (prop.with α b ).

Table S10 :
EGP* candidates proposed during maximization of dipole D with strong constraint on HOMO-LUMO gap ∆ϵ; columns are labeled analogously to TableS7.