Tartarus: A Benchmarking Platform for Realistic And Practical Inverse Molecular Design
Authors: AkshatKumar Nigam
Robert Pollice
Gary Tom
Kjell Jorner
John Willes
Luca Thiede
Anshul Kundaje
Alan Aspuru-Guzik
Venue: NeurIPS 2023 Datasets and Benchmarks Poster (2023)
Keywords: Inverse molecular design
generative models
realistic benchmarking
TL;DR: Realistic benchmarking of generative molecular design models

Abstract: The efficient exploration of chemical space to design molecules with intended properties enables the accelerated discovery of drugs, materials, and catalysts, and is one of the most important outstanding challenges in chemistry. Encouraged by the recent surge in computer power and artificial intelligence development, many algorithms have been developed to tackle this problem. However, despite the emergence of many new approaches in recent years, comparatively little progress has been made in developing realistic benchmarks that reflect the complexity of molecular design for real-world applications. In this work, we develop a set of practical benchmark tasks relying on physical simulation of molecular systems mimicking real-life molecular design problems for materials, drugs, and chemical reactions. Additionally, we demonstrate the utility and ease of use of our new benchmark set by demonstrating how to compare the performance of several well-established families of algorithms. Overall, we believe that our benchmark suite will help move the field towards more realistic molecular design benchmarks, and move the development of inverse molecular design algorithms closer to the practice of designing molecules that solve existing problems in both academia and industry alike.

TARTARUS: A Benchmarking Platform for Realistic And Practical Inverse Molecular Design
AkshatKumar Nigam', Robert Pollice Gary Tom Kjell Jorner John Willes', Luca Thiede", Anshul Kundajc and Alán Aspuru-Guzik 3,4,, Stanford University. Universiity of Groningen. University of Toronto. Vector Institute. Lebovic Fellow. Canadian Institue for Advanced Research. "r. wolliceeruug....o nl
Abstract
The efficient exploration of chemical space to design molecules with intended properties enables the accelerated discovery of drugs, materials, and catalysts, and is one of the most important outstanding challenges in chemistry. Encouraged by the recent surge in computer power and artificial intelligence development, many algorithms have been developed to tackle this problem. However, despite the emer- gence of many new approaches in recent years, comparatively little progress has been made in developing realistic benchmarks that reflect the complexity of molec- this work, we develop set of practical
ular design for real- real-world applications. benchmark tasks relying on physical simulation of molecular systems mimicking real-life molecular design problems for materials, drugs, and chemical reactions. Additionally, we demonstrate the utility and ease of use of our new benchmark set by demonstrating how to compare the performance of several well-established families of algorithms. Surprisingly, we find that model performance can strongly depend on the benchmark domain. We believe that our benchmark suite will help move the field towards more realistic molecular design benchmarks, and move the development of inverse molecular design algorithms closer to designing molecules that solve existing problems in both academia and industry alike.
Introduction
Inverse molecular design, which involves crafting molecules with specific optimal properties M. poses a critical optimization challenge prevalent across various scientific fields. It holds particular significance in chemistry for designing drugs, catalysts, and materials. Each of these scenarios requires a complex balance between multiple desired properties for realistic inverse design. Yet, the majority of comparisons between molecular design algorithms are primarily carried out oversimplified tasks, thereby providing limited insights into their potential performance when confronting intricate design challenges commonplace in chemistry 22.4. While simplistic and cost-effective benchmarks serve as vital tools during the initial stages of algorithm development, facilitating rapid testing and informed design choices, the lack of more realistic benchmarks obstructs their adoption by domain experts in chemistry or materials science. Our work aims to bridge this gap.
In this study, we propose modular suite of design objectives directly inspired by problems regarding new materials, drugs, and chemical reactions. These problems rely on well-established simulation methods in computational chemistry, such as force fields, semi-empirical quantum chemistry, and density functional approximations. We test multiple popular molecular design algorithms on these objectives, offering exemplary performance comparison. The design approaches selected rep- resent some the most important algorithm classes in the field. We also examine the impact of representations the performance of some algorithms. Finally, we provide detailed instructions for
37th Conference on Neural Information Processing Systems (NeurIPS 2023) Track on Datasets and Benchmarks.

---

installation and setup for both the benchmarks and the molecular design algorithms. enabling the cheminformatics and machine learning communities to integrate them into their work and inspire the development of future methods. The core idea is to provide a unified henchmarking framework with a diverse selection of problem domains to assess the generalizability of inverse design algorithms. Notably, the scope of included benchmark domains is to be expanded upon in future work. Altogether, we believe this is a stepping stone towards the widespread adoption of molecular generative models for challenging design problems that permeate the chemical sciences.
2 Background
TARTARUS
Organic photovoltaics
Protein ligands
openbabel
CEP
desigs URROS bis witl sprcle oburnon prourdess
Design sgante for protne degans Souutuo
EdRiE
GDB-13suu
Organic eniiters
Chemical reaction substrates
tbb
Design Erganec receraeiinnso compoN190
Orsign subceibes undegge shubo tydrgges Harcker Explict shuctuo
Dragg bive
crest
10VA0 supp Not ousge
DTP
Datasets
Design benchmarks
Simulation workflows
Cini Venee panclry
Pyscf
SNB-600 60K
aminn
50,000 wyw.syys.soooooo....... werborees
Siigles sijcctee dhuagn loolor weroedda malaaeel repertioes
GoYnnNTee WEREER SEVERINN wobeeryy cIRANNIS
polany
Cesitry
Evaluation Molecular Design Algorithms
Dataset
Properties Ainutaaoon moulss
Molecular Design Algorithms
Simulation workflow
Results Comparison
Seladion Sastapaciic ctitivl
Molecules ProeesGdt condREee
Design benchmark
Set reoaass
Inoogorated peu Noigs dlpellinn Commurison of drinng propored
SUR080O of bolditniai 130h
Figure 1: TARTARUS framework, highlighting its two core elements. Top: Real-world design tasks are defined and paired with simulation workflows and datasets. Bottom: The Evaluation Pipeline: Generative models are conditioned using reference datasets, sampled for structures, and evaluated based on their properties through custom simulation workflows. Model performance is assessed based on the alignment of the acquired properties with predefined objectives.
Generative models play pivotal role in inverse molecular design, yet the progress of benchmark tasks for their performance assessment has been relatively slow [5].. The penalized log P metric, frequently employed in current benchmarks, is characterized by its dependence on molecule size and chain composition, thus limiting its effectiveness 22 [Z-9]. Numerous models have reached saturation point on the task of generating molecules with maximum QED values [10144) The GuacaMol benchmark suite, with its array of goal-oriented objectives (3, often similarly results in near-perfect scores, thus obscuring useful comparisons [15[77]. This concern was traced back to unlimited property evaluations within the suite, with imposed limits revealing significant performance disparities and underscoring the importance of evaluation procedures explicitly accounting for evaluations [18]. The MOSES benchmark suite evaluates the ability of models to propose diverse compounds that mimic the distribution of training datasets (4, However, the emergence of SELFIES and rudimentary algorithms have rendered these tasks trivial 19121]
The recent trend of using molecular docking benchmark for generative models includes tasks like designing molecules complying with Lipinski's Rule of Five (22). However, this often favors highly reactive and unstable molecules 77. Additionally, molecular docking scores have recently become popular metric 223 (2). However, they were only recently incorporated into systematic benchmark platforms. In particular, the Therapeuties Data Commons (TDC) does provide a comprehensive benchmark platform, but crafting sensible, synthesizable, stable structures with low binding affinity estimates remains formidable challenge (28]. The DOCKSTRING package supplies extensive datasets and benchmarks for variety of machine learning strategies in drug design [29]. It constitutes an advance over previous methodologies by offering computational workflows and realistic molecular design benchmarks. However, the breadth of molecular design task extends beyond drug discovery. Accordingly, TARTARUS aspires to cover a more comprehensive set of molecular design problems, to include important structural domains previously neglected.

---

Molecular design algorithms have the potential to accelerate the discovery of materials for the benefit of humankind. However, they can just as well be used design hazardous materials intentionally. While this pertains primarily to the molecular design algorithms themselves, benchmarking can also provide a contribution in that regard by identifying robust and generally applicable tools. However, first, we would like to emphasize that none of the molecular design benchmarks developed have direct link to the development of potentially hazardous materials. Additionally, the molecular design algorithms employed in this work propose merely structures but do not provide any instructions for their synthesis making their malicious use unlikely.
This study introduces four molecular design benchmarks, each with curated dataset of reference molecules with predicted properties, as depicted in Figure [T] The potential availability of reference data from laboratory experiments for at least subsets of our curated datasets was not considered. The goal to present authentic benchmarks that aid discovering structures optimized for specific target properties across diverse applications. To evaluate models with TARTARUS, we recommend training the generative model on the provided dataset, using 80% for training and 20% for hyperparameter op- timization. The model should then propose structures for evaluation by the corresponding benchmark task. Structure optimization is initiated using the best reference molecule from the respective dataset. A constrained approach is followed, limiting the population size, iterations, and total proposed compounds to 5,000, with a capped runtime of 24 hours. Each optimization run is repeated five times for increased robustness and reproducibility. This cespurcce-oostrainnee comparison approach, we argue, is crucial for impartial method comparison and should be a community standard. The parameters and settings used for each model are detailed in the Supporting Information.
3 Results
In this section, we describe the molecular design benchmarks we developed, along with their reference datasets. We also present the results for all the generative models considered, namely: REINVENT (30). SMILES-VAE (31), SELFIES VAE, MoFlow (2), SMILES-LSTM-HC [33 B). SELFIES-LSTM-HE GB-GA 34, and JANUS 77. We would like to emphasize that this is far from comprehensive list of all molecular design algorithms, but rather a small selection of various algorithmic approaches spanning the field. Moreover, the results provided heavily depend on both the model hyperparameters and the benchmark settings, such available computing resources and number of property evaluations. For a comprehensive outline and discussion of the property simulation workflows developed in this work, refer the reader the Additional Results of the Supporting Information.
3.1 Design of Organic Photovoltaics
Euwe
Conformer search
Geometry optimization
Molecule
A.W HOWA LIMI
Scharber model
PCE
Vsita' provesiry form
eptimited geornweyy Piasl
U Alchicabor dipiis evcitet
Calorabbs ottees anagges Dewibe perforsanco
Figure 2: Overview of the property simulation workflow for the design of organic photovoltaics benchmarks. The code accepts SMILES string, generates initial Cartesian coordinates with Open Babel, and performs conformer search and geometry optimization with crest and xtb. Finally, single point calculation at the GFN2-XTB level of theory provides the HOMO and LUMO energies, the HOMO-LUMO gap and the molecular dipole moment. The power conversion efficiency (PCE) is computed from the simulated properties based on the Scharber model.
Our first benchmark objective set comprises two individual tasks inspired by organic photovoltaic (OPV) design (35). The development of organic solar cells (OSCs) garners broad interest due their potential to replace and expand upon the applications of predominant inorganic devices [3 39]. Key properties of photoactive materials for OSCs are their HOMO and LUMO energies (i.e., ΣнOMO and EYUMO. respectively). Importantly, these properties of interest can be directly estimated using standard ground state electronic structure simulation. Thus, we defined two benchmark objectives that are directly inspired by previous high-throughput virtual screening efforts in the literature aimed
3

---

at finding new photoactive materials for OPVs. The first task is based on the Harvard Clean Energy Project (CEP)
40 and corresponds to the design of a small organic donor molecule be used acid methyl ester (PCBM) as acceptor in a bulk heterojunction device The second task based on follow-up work of the CEP and corresponds to the design of a small organic acceptor molecule to be used in bulk heterojunction devices with poly(N- (PCDTBT) as donor [43. In both cases, the objective function is a combination of (1) the estimated power conversion efficiency (PCE) of the corresponding device, which is based on the Scharber model for single junction OSCs 41 42] and derived from the results of electronic structure calculation of the proposed molecule, and (2) the synthethic accessibility of the corresponding structure as quantified by the SAscore
[7
with pheelpheel-l-bbbueee
42 37
40
90-neptadccanyl-2. 1aaprddccaan--2---aaaaaoll---- 7--crbaazlee-llt555---4-4 -,00-,,,,,,,,, 1---44,70--ii22 Aodi-2-thienyl-20. 10..30-enzothhaaiizollll
44
We also provide a reference dataset for training generative models. It is sub- set of the Harvard Clean Energy Project Database (CEPDB) [40), originally en- compassing 2.3 million organic com- pounds. Our subset, CEPsUB. consists of approximately 25,000 molecules. We implemented robust property simula- tion workflow for these benchmark tasks. This workflow accepts the SMILES string of proposed molecules as input [45, 46). generates initial guess of the corre- sponding Cartesian coordinates, performs conformer sampling using the iMTD-GC [8] as implemented in
Table 1: Results for the organic photovoltaics bench- marks. Models are trained subset of the Harvard Clean Energy Project Database. Results are provided as mean and standard deviation of the best target objective values that are obtained in five independent runs in the form mean + standard deviation. The top-performing molecule in the training set is also given ("Dataset"), Met- rics: PCEPCBM: PCBM power conversion efficiency; PC PCEссdтв:: PCDTBT power conversion efficiency; SAscore: synthetic accessibility score.
PCErcus SAwore
PC PCEpory ur SAscuce
Dataset
7.57
31.71
SMILES-VAE SELHES-VAB
7.44 + 0.28
10:23 + 11.14
workflow (7
0.66
4:
crest 49 at the GFN-FF level of the-
MoPlow
7.08 0..1
29.81 0.37
SMILES-LSTM-BC SEI HELMEE--STT--C
6.69 :: 0.40
31.79 1 0.15
ory [50 52 and performs geometry opti-
7 1 0.41
30.71 + 1.20
mization of the lowest energy conformer obtained at the GFN2-XTB level of the- ory [53 54 The predicted properties are taken from a final single point energy cal- culation of the converged geometry. An overview of the simulation workflow and the corresponding benchmark objectives is provided in Figure
7,48 0.11
GB-GA
7.78 = 0.02
30.24 1 0.80
JANUS
7.59 4 0.14
31.34 4 0.74
The results of the two benchmarks are summarized in Table 1 We observe that GB-GA shows best performance on the first benchmark task (7.78), while SMILES-LSTM-HE shows best performance on the second (31.79). While most models are capable to propose molecules with marginally improved PCEs, they are not able to both improve PCEs and diminish SAscores. When looking at the best molecules that the molecular design algorithms generated (cf. Supporting Information), is interesting to see that the various models seem to produce structures with essentially equivalent core motifs. Additionally, it is apparent that, due to the inclusion of the SAscore in the design objectives, the proposed structures are more likely to be stable and synthesizable.
3.2 Design of Organic Emitters
AE(S-TT Singe PAF
Conformer search
Geometry optimization
Molecule
Excited states
Sa(S..-S.) ahargh)
AE(S.SS))
Whhe geamaty Nom
CREST starch
XTB X7B-qpainaee goonsdry
Pina PYSCF singht
Figure 3: Overview of organic emitter design workflow. Starting with a SMILES string, the code executes conformer search and geometry optimization via xtb. TD-DFT single point calculation with pyscf extracts the singlet-triplet gap, oscillator strength, and vertical excitation energy.
The next set of benchmarks is inspired by the design of emissive materials for organic light-emitting diodes (OLEDss), which received significant attention recent years (55 after the discovery of thermally activated delayed fluorescence (TADF) (60]. Their main applications are digital screens

---

and lighting devices [55]. Improving the efficiency of organic emissive materials relying on TADF can be achieved by minimization of their singlet-triplet gaps [6]. Additionally, fluorescence rates need to be increased which corresponds to maximizing the oscillator strength between the first excited singlet and the ground state [[]]. Furthermore, blue OLEDs pose a challenge as they typically suffer from reduced internal quantum efficiencies and device lifetimes [61, (5]. Accordingly, we propose benchmarks targeting the design of emissive materials for OLEDs. The first task is to minimize their singlet-triplet gaps. The second task is to maximize their oscillator strengths. In a third objective, the goal is to combine the first two tasks and, additionally, keep the vertical excitation energies between the ground state and the first excited singlet state in the energy range of blue light. In these tasks, the SAscore of proposed molecules needs to be smaller than or equal to 4.5 for a high fitness.
For training the generative models, we selected a subset of GDB-13 13 [2), which originally consists of more than 970 million organic molecules with up to 13 non-hydrogen atoms, and extracted approximately 380,000 molecules, that consist to a high degree of organic molecules with extended planar conjugated systems of double bonds and lone pairs, i.e., --systems, via a comprehensive set of filters (cf. Supporting Information). The property simulation workflow consists of conformer sampling using the iMTD-GC workflow [47 [S). as implemented in crest [49], at the GFN-FF level of theory 50H52.. followed by geometry optimization of the lowest energy conformer obtained at the GFNO-XTH MEBO-XTBBIO level of theory [53] 63 64) and an excited state single point calculation via TD-DFT at the B3LYP/6-31CH level of theory 70 (Details in the Supporting Information). The entire
simulation workflow is illustrated in Figure
The performance of all the generative models employed on the set of organic emitter design bench- marks is summarized in Table and compared directly to the best baseline fitness values extracted from the dataset. We observe that only JANUS, GB-GA, and the SELFIES- VAE successfully generate compounds that are comparable to or improve upon the best molecules in the training dataset and the lowest singlet-triplet gap with a value of 0.008 eV was found by JANUS. However, this value is only an insignificant improvement over the best reference structure (0.020 eV). These three models also achieved larger oscillator strengths than observed in the training dataset, with GB-GA achieving the highest average oscillator strength with a property value of 2.16. The SELFIES VAE also achieved the highest fitness in the multi-objectiv benchmark task, with an average value of 0.17. All the remaining models were unable to outperform the best molccules in the dataset (-0.04). Finally, the best structures proposed by the generative models in these benchmark tasks (ef. Supporting Information) illustrate that some of them possess reactive structural moieties which is most likely due to the absence of explicit inclusion of stability or synthesizability in the objective functions.
Table 2: Results for the organic emitters design benehmark objectives. Models are trained on a subset of the GDB-13 13 dataset. Results are provided as mean and standard deviation of the best objective values from five independent runs in the form mean + standard deviation. The top-performing molecule in the training dataset is also included ("Dataset"). Metrics: A Ti ): singlet-triplet gap; f12: S1 and So transition oscillator strength; A (F So and Si vertical excitation energy.
AE(S, -Ti- Ti
f13
-V12- AE(S (SE(S-))) -3.
DATASET
0.020
2.97
-0.04
SMILES VAE
0,071 0.003
0.50
+ 0.27
-0.57 + 0.33
SELFIES-VAE
0.016
0.36
0.331 0.31
0.17 + 0.10
MOFLOW
0.001
0.81 0.11
-0.04 + 0.06
-MILES-LSTM-HC
0.015
1,00 0.01
-0.24 + 0.01
FIESLLSMMHH
0.013 0.003
1.00 0.01
-0.24 0.01
REINVENT
0.014 0.003
0.18
-0.15 0.05
GB-GA
0.012 0.002
2.14 0.45
0.07 0.03
JANUS
0.008 0.001
2.07 1 0.16
0.02 1 0.05
3.3 Design of Protein Ligands
We also wanted to include molecular design objectives relevant for medicinal chemistry. In recent years, deep generative models have experienced strong increase in popularity and adoption for drug design as they promise to accelerate discovery campaigns leading to significant cost reductions, and success stories were reported on in the literature 771.74. Therefore, we decided to develop objectives that directly targets the design of ligands for specific proteins based on molecular docking simulations. We set up three benchmark tasks, each aimed at different protein. As targets, we selected ISYH,
5

---

Table 3: Performance metrics for protein- -igann design benchmarks, based on models trained on a subset of the DTP Open Compound Collection. Metries show mean and standard deviation of optimal target objective values over five independent runs (mean + standard deviation). "Dataset" and "Native Docking" values represent the top- perforrrn molecule in the training dataset and the original ligands in their crystal structures, respectively. AEx denotes docking score for protein target X, and SR reflects the success rate for molecules passing structural filters.
ISYH
6Y2F
4LDE
SR
ABsY
AELLDE
SR
QUICK
QUICKVINA?
SMINA
QUICKVINAZ
SMINA
NATIVE DOCKING
-10.2
-10.5
100.0%
-4.9
-5.3
0.0%
-11.6
-12.1
100.0%
DATASET
-10.2
100.0%
-8.2
-8.2
100.0%
-12.2
-13.1
100.0%
SMILES-VAL VAH
0.7
0.6
12.3%
-
-8.9
13.1%
-10.7
+11.1
12.6%
SELFIES VAE
-.3
0.1
34.8%
0.5
-.4
38.9%
-11.1
-12
38.9%
MoFLOW
0.4
0.3
35.5%
9.4
-.3
35.9%
-12.1
0.3
36.2%
SMILES-LSTM-RC
0.6.6:006
-11.1+00
71.7%
0.6.0
-10
-12.4003
0.
73.9%
SELFIES-LSTM-HC
0.5
-12
73.2%
-10.44044
0.6.6
71.5%
-13.7
0.5
75.6%
REINVENT
-11.8 11.81..4
-12.1
77.8%
+..31
1...4.0.3
76.8%
-12.8 0.2
-33.7
4 0.5
76.8%
0.5
+12
72.6%
+.0.9+0..2
73.9%
-3.8
0.4
71.4%
JANUS
-11.220.4
-11
68.4%
-11.3.0..3
-11.9
70.4%
-12.3 0.2
13.6
63.3%
an ionotropic glutamate receptor associated with neurological and psychiatrie diseases such as Alzheimer's, Parkinson's and epilepsy [75] 6Y2F, the main protease of SARS-CoV-2 responsible for translation of the viral RNA of the SARS-CoV-2 virus [76), and 4LDE, the @2-adrenoceptor GPCR receptor spanning the cell membrane and binding adrenaline, hormone which mediates muscle relaxation and bronchodilation [77]. In all cases, crystal structures of the proteins bound to a ligand available in the protein data bank (PDB) [78] and the objective is to minimize the docking score of a proposed molecule to one these protein targets.
For training, we started from the De- velopmental Therapeutics Program (DTP) Open Compound Collection 80].
Molecular Docking
set of about 250,000 molecules tested for treatment against cancer and the acquired immunodefi- ciency syndrome (AIDS) (81), and se- lected structures passing the struc- ture constraints leading to around 150,000 structures referred to as DATASET. Molecules included in this collection can be ordered for free, except for shipping costs, from the DTP 82). The simulation workflow for these benchmark tasks is initiated with the SMILES string of the pro- posed molecule and creates an initial guess of the corresponding Cartesian coordinates with openbabel. Sub- sequently, we use QuickVing [83) introduce the newly proposed com- pound in the protein binding pocket, and perform a docking pose search. Following recent literature recommen- dations 84 85), the top compounds from each run were re-scored for improved accuracy using
79
Molecule
A Docking score target
Initial geometry from openbabel
QuickVina2 pase sampling Smina re-scoring
docking site
Protein targets
GYZF
1SYH
4LDE
Figure 4: Overview of the computational workflow for pro- tein ligands design benchmarks. The process accepts SMILES string for protein targets ISYH, 6Y2F, and 4LDE, generates Cartesian coordinates using Open Babel, and sam- ples docking poses via smina. The lowest docking score is returned the target property.
smina 86 (see Figure
The results are compiled in Table which includes the averages of the top docking scores achieved by each model. We also present the highest docking scores from the reference dataset for each of the tasks (DATASET Table along with the docking scores the original ligands bound to the three target proteins in their respective crystal structures ("Native Ligand" in Table Moreover, Table 3 indicates the percentage of sampled molecules that pass standard structure filters employed drug discovery (see SR in Table 3). Notably, the native ligand for 6Y2F fails to pass these constraints, illustrating that these constraints only evaluate the likelihood of a structure being both stable and synthesizable. Across all protein targets,SMILES. EAEILLEES--EEEEEEEE struggles to suggest structures with better

---

docking scores than in the reference dataset (-10.0/-10.4 for ISYH, -8.3/-8.9 for 6Y2F, and 10.77-11.1 for 4LDE), and it also yields the lowest success rates (12.3% for 1SYH, 13.1% for 6Y2F, and 12.6% for 4LDE). The SELFIES- VAE appears to perform slightly better for both docking scores (-10.4/-10.9 for ISYH, -9.6/-10.1 for 6Y2F, and 11.1/-11.9 for 4LDE) and success rates (34.8% for ISYH, 38.9% for 6Y2F, and 38.9% for 4LDE), but not significantly so. No single model consistently achieves the highest docking scores across all three targets. REINVENT attains the lowest average docking score for ISYH (-11.8/-12.1) JANUS for 6Y2F (-11.3/-11.9) and GB-GA for 4LDE (-12.9/-13.8) Finally, due to the integration of our filters, the best designs in this set benchmarks largely correspond to stable and potentially synthesizable compounds.
3.4 Design of Chemical Reaction Substrates
Reactant
APP Frs En
Crarg
SEAM optimization
Transition state
Nedeers senkiie calonarr
theer colivealee pronnete)
VIGAOW
The www.eoooiooooooo GETTIGATT
Product
DE: Ep-
En
Sraige pickuct
Figure 5: Overview of the workflow for designing chemical reaction substrates. Starting with SMILES string, the code optimizes reactants and products, generates guessed transition state via SEAM optimization, and conducts constrained conformational sampling. Reaction energy and approximate SEAM activation energy are then extracted.
Table 4: Results for chemical reaction substrate design benchmarks. Models, trained benchmark dataset generated with STONED-SELFIES cycles, yield mean and standard deviation of optimal objective values over five runs (mean + standard deviation). Baselines include the best molecule in the training dataset ("Dataset") and the parent unsubstituted substrate ("Parent Substrate"). Metrics: A F Activation energy of the reaction; AF,: Reaction energy.
A
AET AF,
AE AE,
PARENT SUBSTRATE
83.16
0,00
$5.16
-85.16
DATASET
64.94
-34.39
56.48
-95.25
SMILES-VAE
76.81
0.25
-17.96
71.01
0.62
-90.94
FIES-VAE
72.45 = 3.79
-10.45
72.05 0 0.00
.87.82 2.13 2.13
MoFLOW
63.21
-92.82
MMILES-LSTM-HE
59.64 = 4.10
-31.03
16.15
71.81 = 1.56
-91.58 2.14
63.17 4.34
-21.02
4.95
68.06 5.74
.96.59 4.59
REINVENT
68.38 2,00
24.35
6.46
55.25
5.88
.94.52 1 1.20
GR-GA
56.04 3.07
.41.99
5.76
45.20
6.78
.100007
1.38
47.56
445.37
7.90
Developing new chemical reactions and finding new catalysts for existing ones are important goals drive innovations drug and material discovery [87]. and move towards more sustainable pro- duction 88 Due the importance of chemical reactions and the lack of reliable molecular design benchmarks related to reactivity, we looked for alternative approaches allowing us to model TSs in reliable way with increased efficiency. We decided to use the SEAM force field method, which combines the force fields two minima directly connected via the intrinsic reaction coordinate to the TS of interest generating effective TS force field [89].
With this method, chose model the intramolecular concerted double hydrogen transfer reaction of seseesuuuoorrosssnss which is a dyotropic reaction of type II [90]. We use this reaction to define an inverse molecular design benchmark that targets modification of the substrate and product structures to alter reactivity. As main target properties defining reactivity in this system, we selected corresponding activation and reaction energies. The first objective to minimize the activation energy and corresponds to making the reaction faster, a common target when developing catalysts or reagents. Our second objective focuses on the thermodynamic driving force and aims to minimize the reaction energy making the respective product maximally thermodynamicall favorable. Both the

---

third and the fourth task combine two properties into one objective function. The third corresponds to finding substrates that lead to both a fast and thermodynamically favorable reaction. The fourth is about finding substrates that cause both a slow and a thermodynamically favorable reaction.
Like for the other benchmarks, we provide a reference dataset that be used for training the generative models. As this benchmark requires molecules that contain the syn=sesquinorborrene motif, we needed create the new SNB-60K dataset with approximately 60,000 molecules. The simulation workflow starts with a check of the required structural constraints (vide supra) in the proposed substrate SMILES string. If any constraint is violated, the workflow will be aborted and an extremely bad objective value of -0,000 returned. If all constraints are satisfied, initial guess Cartesian coordinates will be generated. Subsequently, combination of conformer searches and constrained optimizations generates the reactant and product structures required to set up the SEAM force field and perform TS optimization. Finally, conformer search of the TS is carried out, followed by subsequent SEAM force field optimization, and the reaction and activation energies are derived from the lowest energy structures of reactant, TS and product. The workflow is shown in Figure [5]
The performance of the generative models considered on the molecular reactivity benchmarks is provided in Table As similarly seen on the protein ligand design benchmarks, we observed that both the VAE models failed to surpass the best values in the dataset (e.g.. SMILES-VAB 76.81 versus 64.94 in the dataset for AE', -10.96 versus -34.49 in the dataset for A E,. 71.01 versus 56.48 in the dataset for AE + AE, and -90.94 versus -95.25 in the dataset for E + A, Compared to the VAEs, both the LSTM-HC models, MoFlow, and REINVENT perform better but they are not able to improve upon the best molecules from the dataset for all four benchmark objectives. In particular, SELFIES-LSTMMHE only reaches -21.00 21.02 versus -3..49 the dataset for A Fr, and 68.06 versus 56.48 in the dataset for A + A,, Only JANUS and GB-GA consistently propose structures that outperform the best reference compounds. The best structures proposed in these benchmark tasks (cf. Supporting Information) generally show that the majority of structures is reasonable in terms of stability. In addition, the inclusion of SAscore constraints in the two combined objectives seems to particularly result in smaller molecules as molecular size tends to correlate with the SAscore.
3.5 Model Timing Comparison
Finally, one overlooked consideration for choosing generative models is the computation time needed for pre-conditioning the model based a reference dataset and for proposing new candidates. The main reason we conducted this benchmark is our belief that shorter pre-conditioning and structure generation translates increased usage of molecular design algorithms. While it is a minor overhead in the case of time-limiting property evaluation, it can constitute a significant entrance barrier for testing a model in benchmark tasks with comparably fast evaluations before using it on the actual problem of interest. We used three different metrics for comparison based on the reference dataset of the design of protein ligands. First, we evaluated the time the models need for pre-conditioning (or training) when provided with the corresponding reference dataset. For that, the models have access to 24 CPU cores and GPU. Second, as access to GPUs can sometimes be limited, we wanted to evaluate the impact of GPU usage on the model pre-conditioning time. However, for some of the models tested, training times were too prohibitive for benchmarking when no GPU was provided. Thus, to measure the impact of using a GPU for training, instead of the full time, we evaluated single epoch times both with and without access to a GPU. Third, we determined the time required by the algorithms to propose 10,000 unique structures. All timing results are illustrated in Figure 6
Overall, we find that both VAEs needed by far the longest training time, both taking considerbbyy longer than 9 hours, with the SEI FIELFIES-VE VAH requiring somewhat more time than the SMILES- VAÉ. In contrast, both LSTM-HC models only need approximately 2 hours with no statistically significant difference between using SMILES and SELFIES. Compared to the other deep generative models, both MoFlow and REINVENT show by far the fastest training times finishing already after less than 1 hour in our benchmark. When looking at the single epoch times with and without GPU usage, we find that the relative order with respect to timing is comparable regardless of whether GPU is utilized. It is particularly notable that both the VAEs have single epoch times of approximately 8 hours which leads to an estimated training time without GPU of about 25 days, which we decided to be too prohibitive to perform full training. All other deep generative models have single epoch times significantly faster than 1 hour leading to estimated training times below day. The fastest models in that respect are again MoFlow and REINVENT with single epoch times of only several
8

---

trained minutes. Notably, MoFlow has a faster epoch time than REINVENT with GPU usage, but slower one without. Nevertheless, the total pre-conditioning times of both GAs, which are provided in the diagram depicting single epoch times (cf. Figure (C), is still shorter which only surmounts to several seconds without the use of a GPU.
4 Conclusion and Outlook
We developed TARTARUS, a modular set of realistic molec- ular design objectives represen- tative of common problems in the domains of organic materi- als, drugs, and chemical reac- tions to be used as general goal- directed benchmarks for gener- ative models. Most tasks are directly inspired by design ob- jectives pursued in the litera- ture. Hence, molecules have potential value in real applications. We observed that none of the representative generative models selected con- sistently outperformed the oth- ers across all the design tasks. This is due a seemingly strong dependence of the performance of some models on the specific benchmark domain and suggests the importance of conducting di- verse molccular design bench- marks that reflect actual prob- lems pursued in chemistry to evaluate the generalizability of inverse molecular design algorithms across domains. Furthermore, it shows that current molecular design algorithms have significant room for improvement in terms generality, particularly deep generative models. Thus, more sophisticated approaches do not guarantee stronger inverse design performance, especially when not tested on the problem domain of interest.
A)
Training Time
B) 10k Motecule Sample Time
SE. FIES LT HC:
SWI ES STM
Time (secondss
SNILES CEILES-LSTHIE LSTH-K
Time (secondis)
SWILESS-V
Stt Pi-------
REDWENT
RES RESWENT
Moflow
Epoch Time (CPU)
D)
Time (GPU)
STM--T
Time (seconds)
REINIENT
Time (sscands))
SMILES LST HC
SWILESS VAE
SE.PIES
5 SMILLS
GA
SELFIES
MoToow
SMSLES
JNNUS
Figure 6: Results of model timing comparisons. Models are subset of the DTP Open Compound Collection. The benchmark metrics consist of model training time (or pre- condition time) with the of GPU (A), single epoch times during training both with (D) and without a GPU (C), and the time required to sample 10,000 structures (B). Results are provided as mean (main bar) and standard deviation (error bars) of the metrics obtained in five independent runs.
We note that relatively simple approaches like genetic algorithms that do not need large reference datasets show more consistently good performance across benchmarks, whereas the VAE approaches, that strongly rely on the available training data, show limited performance across multiple bench- marks. However, more comprehensive investigations are required confirm this hypothesis and identify the origins of the limited performance of some of the approaches tested. Moreover, we believe that demonstrating strong performance of computer algorithms applied realistic design objectives will inspire researchers across fields to adopt these techniques in their workflows eventually leading to mainstream use. Accordingly, we believe that TARTARUS will act as a stepping stone for advancing deep generative models and for inspiring the development of realistic molecular design benchmarks. Notably, we would like to emphasize that TARTARUS does not replace the development of more specialized molecular design benchmarks that are tailored to specific domains of interest and necessarily follow a different design workflow, particularly in drug design [9I-94). We believe that such specialized benchmarks complementary TARTARUS and are equally important to assess the practicality and performance of molecular generative models.
Nevertheless, this first work on TARTARUS still leaves room for further developments. As is common real-world molecular design, objectives need to be revised when undesired structures are promoted and new property requirements can emerge. In particular, the relative importance of multiple target properties sometimes needs to be adjusted and the consideration of additional properties becomes necessary. Furthermore, the benchmark domains currently covered are far from comprehensive,

---

providing opportunities for the development of additional design tasks in the near future. Deep generative models capable design 3D molecular structures directly are currently not well supported. as proposed conformers would be ignored. Extending TARTARUS in that regard requires including geometries our reference datasets. Both for the protein ligand and the chemical reaction substrate design benchmarks, specialized geometries are needed. Hence, it would not be enough to train 3D) generative models suggesting ground-state gas-phase geometries. The chemical reaction substrate design benchmarks require multiple geometries as both two ground state geometries and one transition state geometry are simulated. Consequently, most current 3D generative models cannot fulfill this requirement out of the box and could thus not replace established global geometry optimization strategies. Finally, the results of 3D generative models require dedicated performance comparisons molecular properties are intrinsically tied to their 3D structures, For meaningful comparison to the other molecular design algorithms considered, it is necessary to separate 3D structure generation from molecular connectivity generation. This requires introducing separate conformer generation benchmarks for each of the benchmark tasks. Accordingly, these extensions of TARTARUS will be the subject of future work. Finally, aim launch open molecular design challenges making use of the benchmark codes included in TARTARUS. This will allow us to provide regular measurements of the state of the art in molecular generative models and promote further advances in the field.
regular basis
Data Availability
Code and datasets for running all the TARTARUS benchmarks are provided in our public GitHub repository: hwwp:///ww.....oaee////////////e
(DOI:
For further discussions and collaboration, we invite readers to become part of our Discord community;
Acknowledgements
AK.N. acknowledges funding from the Bio-X Stanford Interdisciplinary Graduate Fellowship (SGIF). R.P. acknowledges funding through Postdoc. Postdoc.Mobility fellowship by the Swiss National Science Foundation (SNSF, Project No. 191127). G.T. acknowledges the Natural Sciences and Engineering Research Council of Canada for support through the Postgraduate Scholarships-Doctore (PGS-D) pro- gram. K.J. acknowledges funding through an International Postdoc grant from the Swedish Research Council (No. 2020-00314). A.A.-G. thanks Anders G. Frøseth for his generous support. A.A.-G, ac- knowledges the generous support of Natural Resources Canada and the Canada 150 Research Chairs program. Some computations were performed on the Béluga and Narval supercomputers situated the Éccole de technologie supérieure in Montreal. We also thank the SciNet HPC Consortium for support regarding the use of the Niagara supercomputer. SciNet funded by the Canada Foundation for Innovation, the Government of Ontario, Ontario Research Fund Research Excellence, and the University of Toronto. Computations were also performed on the Cedar supercomputer situated at the Simon Fraser University Burnaby. In addition, we acknowledge support provided by Compute Ontario and Compute Canada. Some of the computing for this project was also performed on the Sher- lock cluster. We would like thank Stanford University and the Stanford Research Computing Center for providing computational resources and support that contributed to these research results. We thank the U.S. Department of Energy (DOE)/REE//ALLIANN Spectra (AM1.5G) publicly available at https///www...rrroooorroddddeeee
for making the Reference Air Mass 1.5 1.5.html
10

---

References
[1] Benjamin Sanchez-Lengeling and Alán Aspuru-Guzik. Inverse molecular design using machine learning: Generative models for matter engineering. Science, 361(64000:360-365, 2018.
[2] Rafacl Gómez-Bombarelli, Jenmifer N. Wei, David Duvenaud, José Miguel Hernández- Lobato, Benjamín Sánchez-Lengeling Dennis Sheberla, Jorge gguilera-Iparraguirs Timothy D. Hirzel, Ryan P. Adams, and Alán Aspuru-Guzik. Automatic chemical design using a data- driven continuous representation of molecules. arXiv, page 1610.02415v1, 2016.
[3] Nathan Brown, Marco Fiscato, Marwin HS Segler, and Alain C Vaucher. Guacamol: bench- marking models for de novo molecular design. Journal of chemical information and modeling. 19(3):106-1108) 2019.
[4] Daniil Polykovskiy, Alexander Zhebrak, Benjamin Sanchez-Lengeling, Sergey Golovanov, Oktai Tatanov, Stanislav Belyaev, Rauf Kurbanov, Aleksey Artamonov, Vladimir Aladinskiy, Mark Veselov, Artur Kadurin, Simon Johansson, Hongming Chen, Sergey Nikolenko, Alán Aspuru-Guzik, and Alex Zhavoronkov. Molecular sets (moses): A benchmarking platform for molecular generation models. Frontiers in Pharmacology, 11, 2020.
[5] Gisbert Schneider and Uli Fechner. Computer-based de novo design of drug-like molecules. Nature Reviews Drug Discovery, 4(8):649 663, August 2005.
[6] Tony Slater and Dave Timms. Meeting on binding sites: Characterizing and satisfying steric and chemical restraints. University of York, 28-30 March 1993. Journal of Molecular Graphics, ((4):248-251 December 1993.
[7] AkshatKumar Nigan, Robert Pollice, and Alán Aspuru-Guzik. Parallel tempered genetic algorithm guided by deep neural networks for inverse molecular design. Digital Discovery, page Advance Article, 2022.
[8] AkshatKumar Nigam, Pascal Friederich. Mario Krenn, and Alan Aspuru-Guzik. Augmenting genetic algorithms with deep neural networks for exploring the chemical space. In International Conference on Learning Representations, 2020.
[9] Xiufeng Yang, Jinzhe Zhang, Kazuki Yoshizoe, Kei Terayama, and Koji Tsuda. Chemts: an efficient python library for de novo molecular generation. Science and technology of advanced materials, (11)9972-976, 2017.
[10] Jiaxuan You, Bowen Liu, Zhitao Ying, Vijay Pande, and Jure Leskovee. Graph convolu- tional policy network for goal-directed molecular graph generation. In Advances in Neural Information Processing Systems, pages 6410-6421, 2018
[1I] Mariya Popova, Mykhailo Shvets, Junier Oliva, and Olexandr Isayev. Molecularrnn: Generat- ing realistic molecular graphs with optimized properties. arXiv preprint 1905v:1905.13377 2019.
[12] Zhenpeng Zhou, Steven Kearnes, Li, Richard N Zare, and Patrick Riley. Optimization of molecules via deep reinforeement learning. Scientific reports, 9()::-100, 2019.
[13] Chence Shi*, Minkai Xu*, Zhaocheng Zhu, Weinan Zhang, Ming Zhang, and Jian Tang. Graphaf: a flow-based autoregressive model for molecular graph generation. In International Conference on Learning Representations, 2020.
[14] Wengong Jin, Regina Barzilay, and Tommi Jaakkola. Junction tree variational autoencoder for molecular graph generation. In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 2323-2332. PMLR, 10-15 Jul 2018.
[15] Sungsoo Ahn, Junsu Kim, Hankook Lee, and Jinwoo Shin. Guiding deep molecular op- timization with genetic exploration. Advances in neural information processing systems, 13:12008-1202 2020.

---

[16] Jules Leguy. Thomas Cauchy, Marta Glavatskikh, Béatrice Duval, and Benoit Da Mota. Evomol: a flexible and interpretable evolutionary algorithm for unbiased de novo molecular generation. Journal of Cheminformatics, 12(1):1-19, 2020.
[17] Yongbeom Kwon and Juyong Lee. Molfnnde: an evolutionary algorithm for the global optimization of molecular properties and the extensive exploration of chemical space using smiles. Journal of cheminformaties, 13(1):1-14, 2021.
[18] Wenhao Gao, Tianfan Fu, Jimeng Sun, and Connor W Coley. Sample efficiency matters: A benchmark for practical molecular optimization. arXiv preprint arXiv:2206. 124/1, 2022.
[19] Mario Krenn, Florian Häse, AkshatKumar Nigam, Pascal Friederich, and Alan Aspuru- Guzik, Seff-referencing embedded strings (selfies): A 100% robust molecular string representation. Machine Learning: Science and Technology. 1(4):045024, 2020.
[20] Mario Krenn, Qianxiang Ai, Senja Barthel, Nessa Carson, Angelo Frei, Nathan Frey, Pascal Friederich, Théophile Gaudin, Alberto Alexander Gayle, Kevin Maik Jablonka, et al. Selfies and the future of molecular string representations. arXiv preprint arXiv:2204.00055, 2022.
[21] Philipp Renz, Dries Van Rompaey, Jörg Kurt Wegner, Sepp Hochreiter, and Günter Klam- bauer. On failure modes in molecule generation and optimization. Drug Discovery Today: Technologies, 32:55-65 63, 2019.
[22] Tobiasz Cieplinski, Tomasz Danel, Sabina Podlewska, and Stanislaw Jastrzebski. Generative models should at least be able to design molecules that dock well: A new benchmark. Journal of Chemical Information and Modeling, 63(11):3238-322 3247, 2023. PMID: 37224003.
[23] Gisbert Schneider and David E. Clark. Automated de novo drug design: Are we nearly there yet? Angew. Chem. Int. Ed., (8(32):10722-10803 2019.
[24] Woosung Jeon and Dongsup Kim. Autonomous molecule generation using reinforcement learning and docking to develop potential novel inhibitors. Scientific reports, 10(1):22104, 2020.
[25] Sarvesh Mehta, Siddhartha Laghuvarapu, Yashaswi Pathak, Aaftaab Sethi, Mallika Alvala, and U. Deva Priyakumar. Memes: Machine learning framework for enhanced molecular screening. Chem. Sci., 12:11700-1721 2021.
[26] Morgan Thomas, Robert T Smith, Noel M O'Boyle, Chris de Graaf, and Andreas Bender. Comparison of structure-and ligand-based scoring functions for deep generative models: a gper case study. Journal of Cheminformatics, 13(1):1-20, 2021.
[27] Tiago Sousa, João Correia, Vitor Pereira, and Miguel Rocha. Generative deep learning for targeted compound design. Journal of Chemical Information and Modeling, 61(11):5343-531, 2021. PMID: 34699719.
[28] Kexin Huang, Tianfan Fu, Wenhao Gao, Yue Zhao, Yusuf H Roohani, Jure Leskovee, Connor W. Coley, Cao Xiao, Jimeng Sun, and Marinka Zitnik. Therapeutics data commons: Machine learning datasets and tasks for drug discovery and development. In Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 1), 2021.
[29] Miguel Carcia-Ortegge Gregor N. C. Simm, Austin J. Tripp. José Miguel Hernández- Lobato, Andreas Bender, and Sergio Bacallado. Dockstring: Easy molecular docking yields better benchmarks for ligand design. Journal of Chemical Information and Modeling, 0(0),uull, 2022.
[30] Mareus Olivecrona, Thomas Blaschke, Ola Engkvist, and Hongming Chen. Molecular dovoovo design through deep reinforcement learning. Journal of cheminformatics, 911):1-14, 2017.
2017
[31] Rafael Gómez-Bombarelli, Jennifer N. Wei, David Duvenaud, José Miguel Hernández-Lobato Benjamín Sánchez-Lengeling, Dennis Sheberla, Jorge eguilera-lparraguirre Timothy D. Hirzel, Ryan P. Adams, and Alán Aspuru-Guzik. Automatic chemical design using data- driven continuous representation of molecules. ACS Central Science, 4(2):268-277 Jan
2018.
12

---

[32] Chengxi Zang and Fei Wang. Moflow: an invertible flow model for generating molecular graphs. Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining, pages 617-626, 2020.
[33] Marwin HS Segler, Thierry Kogej, Christian Tyrchan, and Mark P Waller. Generating focused molecule libraries for drug discovery with recurrent neural networks. ACS central science, 4(1):120-131, 2018.
[34] Jan H Jensen. A graph- genetic algorithm and generative model/monte carlo tree search for the exploration of chemical space. Chemical science, 10(12):35773377, 2019.
[35] Alexander W. Hains, Ziqi Liang, Michael A. Woodhouse, and Brian A. Gregg. Molecular semiconductors in organic photovoltaic cells. Chemical Reviews, 110(11)66699-635 2010.
PMID: 20184362.
[36] Sean E Shaheen, David S Ginley, and Ghassan E Jabbour. Organic-based photovoltaies: toward low-cost power generation. MRS bulletin, 30(1):10-19 2005.
[37] Johannes Hachmann, Roberto Olivares- Sule Atahan-Evrente Carlos Amador-Bedolle Roel S. Sánchez-Carrera, Aryeh Gold-Parker, Leslie Vogt, Anna M. Brockway, and Alán Aspuri-Guzik. The harvard clean energy project: Large-scale computational screening and design of organic photovoltaics on the world community grid. The Journal of Physical Chemistry Letters, 2(17):2241-225) 2011.
[38] Hongseok Youn, Hui Joon Park, and L. Jay Guo. Organic photovoltaic cells: From performance improvement to manufacturing processes. Small, 11(19):2228-2245 2015.
[39] Moritz Riede, Donato Spoltore, and Karl Leo. Organic solar cells cells-the the path to commercial success, Advanced Energy Materials, 11(1)22002653, 2021.
[40] Johannes Hachmann, Roberto Olivares- Adrian Jinich, Anthony L. Appleton, Mar- tin A. Glood-Forsytte László R. Seress, Carolina SalgndSal Kai Trepte, Sule Atahan- Evrenk, Süleyman Er, Supriya Shrestha, Rajib Mondal, Anatoliy Sokolov, Zhenan Bao, and Alán Aspuru-Guzik. Lead candidates for throughput quantum chemistry the harvard clean energy project. Energy Environ. Sci., 2014.
organic photovoltaics from high-
[41] M. C. Scharber, D. Mühlbacher, M. Koppe, P. Denk, C. Waldauf, J. Heeger, and C. Brabec. Design rules for donors bulk-heterojunction solar towards 10 % energy-conversion efficiency. Advanced Materials, 18(6):789-794, 2006.
[42] Tayebeh Ameri, Gilles Dennler, Christoph Lungenschmied, and Christoph J. Brabec. Organic tandem solar cells: A review. Energy Environ. Sci., 2:347-363, 2009.
[43] Steven A. Lopez, Benjamin Sanchez-Lengeling, Julio de Goes Soares, and Alán Aspuru-Guzik. Design principles and top non-fullerene acceptor candidates for organic photovoltaics. Joule, 1(4):857-870, 2017.
[44] Peter Ertl and Ansgar Schuffenhauer. Estimation of synthetic accessibility score of drug-like molecules based on molecular complexity and fragment contributions. Journal of cheminfor- maties, 1(1):8, 2009.
[45] David Weininger. Smiles, a chemical language and information system. 1. introduction to methodology and encoding rules. Journal of chemical information and computer sciences, 28(1):31-36, 1988.
[46] David Weininger, Arthur Weininger, and Joseph L.. Weininger. Smiles. 2. algorithm for generation of unique smiles notation. Journal of Chemical Information and Computer Sciences, 19 1989.
[47] Stefan Grimme. Exploration of chemical compound, conformer, and reaction space with meta-dynamics simulations based on tight-binding quantum chemical calculations. Journal of Chemical Theory and Computation, 15(5):2847-2862. 2019. PMID: 30943025.
13

---

[48] Philipp Pracht, Fabian Bohle, and Stefan Grimme. Automated exploration of the low-energy chemical space with fast quantum chemical methods. Phys. Chem. Chem. Phys., 22:7169-7199 2020.
[49] Philiip Pracht and Stefan Grimme. ConformerRRotmee Ensemble Sampling Tool. http://gihub.coo//iinnellbllaoi
May 2020.
[50] Sebastian Spicher and Stefan Grimme. Efficient computation of free energy contributions for association reactions of large molecules. The Journal of Physical Chemistry Letters, 11(16):6606-66 11176666666111 2020. PMID: 32787231.
[51] Philipp Pracht, David F. Grant, and Stefan Grimme. Comprehensive assessment of gfn tight- binding and composite density functional theory methods for calculating gas-phase infrared spectra. Journal of Chemical Theory and Computation. 16(11):7044-706) 2020. PMID: 33054183.
[52] Sehastian Spicher and Stefan Grimme. Robust atomistic modeling of materials, organometallic, and biochemical systems. Angewandre Chemie International Edition, 59(36);15665-1567, 2020.
[53] Stefan Grimme, Christoph Bannwarth, and Philip Shushkov. A robust and accurate tight- binding quantum chemical method for structures, vibrational frequencies, and noncovalent interactions of large molecular systems parametrized for all spd-block elements (z 1-86). Journal of Chemical Theory and Computation, 13(5):199-20099 2017. PMID: 28418654.
[54] Christoph Bannwarth, Sebastian Ehlert, and Stefan Grimme. GFN2-xTB An accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependen dispersion contributions. Journal of Chemical Theory and Computation, 15(3):1552-1671 March 2019.
[55] Michael Y. Wong and Eli Wesmann-Clman Purely organic thermally activated delayed fluorescence materials for organic light-emitting diodes. Advanced Materials, (9teriaall,9222)116666 29(22):1605444, 2017.
[56] Zhiyong Yang, Zhu Mao, Zongliang Xie, Yi Zhang, Siwei Liu, Juan Zhao, Jiarui Xu, Zhen- guo Chi, and Matthew P. Aldred. Recent advances in organic thermally activated delayed fluorescence materials. Chem. Soc. Rev, 46:915-1016 2017.
[57] Yuchao Liu, Chensen Li, Zhongjie Ren, Shouke Yan, and Martin R Bryce. All-organic thermally activated delayed fluorescence materials for organic light-emitting diodes. Nature Reviews Materials, 3(4):1-20, 2018.
[58] AkshatKumar Nigam, Robert Pollice, Pascal Friederich, and Alán Aspuru-Guzai Guzik. Artificial de- sign of organic emitters via genetic algorithm enhanced by deep neural network. ChemRxiv, 2023.
[59] Robert Pollice, Benjamin Ding, and Alán Aspuru-Guzik. Rational design of organic molecules with inverted gaps between first excited singlet and triplet. ChemRxiv, 2023.
[60] Hiroki Uoyama, Kenichi Goushi, Katsuyuki Shizu, Hiroko Nomura, and Chihaya Adachi. Highly efficient organic light-emitting diodes from delayed fluorescence. Nature, 492(7428):234 2012.
[61] Rafael (ö6mez-Bombarde Jorge Aguilera-Iparraguires Timothy Hirzel, David Duvenaud, Dougal Maclaurin, Martin A Blood-Forsythe, Hyun Sik Chae, Markus Einzinger, Dong-Gwang Ha, Tony Wu, et al. Design of efficient molecular organic light-emitting diodes by a high- throughput virtual screening and experimental approach. Nature materials, 15(10): 120-1 2016.
[62] Lorenz C Blum and Jean-Louis Reymond. 970 million druglike small molecules for virtual screening in the chemical universe database gdb-13, Journal of the American Chemical Society, 131(25):8732-877) 8733, 2009.
14

---

[63] Philipp Pracht, Eike Caldeweyher, Sebastian Ehlert, and Stefan Grimme. A Robust Non-Sell Consistent Light-Binding Quantum Chemistry Method for large Molecules. ChemRxiv, June 2019.
[64] Christoph Bannwarth, Eike Caldeweyher, Sebastian Ehlert, Andreas Hansen, Philipp Pracht, Jakob Seibert, Sebastian Spicher, and Stefan Grimme. Extended tight-binding quantum chemistry methods. WIREs Computational Molecular Science, 11(2):e 1493, 2021.
[65] A. D. Becke. funntty-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A, 18:3098-31000 Sep 1988.
[66] Chengteh Lee, Weitao Yang, and Robert G. Parr. Development of the colle-salvetti correlation- energy formula into a functional of the electron density. Phys. Rev. B, 37:785-789, Jan 1988.
[67] Axel D. Becke. Density-functione thermochemistry. iii. the role of exact exchange. The Journal of Chemical Physics, 98(7):5648 56522, 1993.
[68] R. Ditchfield, W. J. Hehre, and A. Pople. Self-consistent molecular-orbital methods. ix. an extended gaussian-typ basis for molecular-orbital studies of organic molecules. The Journal of Chemical Physics, 54(2):724-728, 1971.
[69] W. J. Hehre, R. Ditchfield, and J. A. Pople. Self-consistent molecular orbital methods. xii. further extensions of gaussian-type basis sets for use in molecular orbital studies of organic molecules. The Journal of Chemical Physics, 5665)):2577-2261 1972.
[70] Praveen C Hariharan and John A Pople. The influence of polarization functions on molecular orbital hydrogenation energies. Theoretica chimica acta, 28(3):213-222, 1973.
[71] Alex Zhavoronkov, Yan A Ivanenkov, Alex Aliper, Mark S Veselov, Vladimir A Aladinskiy. Anastasiya V Aladinskaya, Victor A Terentiev, Daniil A Polykovskiy, Maksim D Kuznetsov, Arip Asadulaev, et al. Deep learning enables rapid identification of potent ddrl kinase inhibitors. Nature biotechnology, 37(9):1038 1040, 2019.
[72] Kit-Kay Mak, Madhu Katyayani Balijepalli, and Mallikarjuna Rao Pichika. Success stories of ai in drug discovery where do things stand? Expert Opinion on Drug Discovery, 17(1):79-92. 2022. PMID: 34553659.
[73] Atanas Patronov, Kostas Papadopoulos, and Ola Engkvist. Has Artificial Intelligence Impacted Drug Discovery?, pages 153-176. Springer US, New York, NY, 2022.
[74] Yunan Luo, Jian Peng, and Jianzhu Ma. Next decade's ' ai-based drug development features tight integration of data and computation. Health Data Science, 2022, 2022.
[75] Anne Frandsen, Darryl S Pickering, Bente Vestergaard, Christina Kasper, Bettina Bryde Nielsen, Jeremy R Greenwood, Giuseppe Campiani, Caterina Fattorusso, Michael Gajhede, Arne Schoushoe, et al. Tyr702 an important determinant agonist binding and domain closure of the iigand-bindigg core of glur2. Molecular pharmacology. 67(3):703-713, 2005.
[76] Linlin Zhang, Daizong Lin, Xinyuanyuan Sun, Ute Curth, Christian Drosten, Lucie Sauer- hering, Stephan Becker, Katharina Rox, and Rolf Hilgenfeld. Crystal structure of sars-cov-2 main protease provides a basis for design of improved a-ketoamide inhibitors. Science, 368(6489):409444 2020.
[77] Aaron M Ring, Aashish Manglik, Andrew C Kruse, Michael D Enos, William I Weis, K Christo- pher Garcia, and Brian K Kobilka. Adrenaline-activaaeed structure of 3 2-adrenoceptor stabi- lized by an engineered nanobody. Nature, 502((742)):755557777 2013.
[78] Helen M. Berman, John Westbrook, Zukang Feng, Gary Gilliland, T. N. Bhat, Helge Weissig, Ilya N. Shindyalov, and Philip E. Bourne. The Protein Data Bank. Nucleic Acids Research, 28(1):235 245,,222,001 01 2000.
[79] Johannes H. Voigt, Bruno Bienfait, Shaomeng Wang, and Mare C. Nicklaus. Comparison of the nci open database with seven large chemical structural databases. Journal of Chemical Information and Computer Sciences, 41(3):702-712, 2001.
15

---

[80] Wolf-Dietries Ihlenfeldt, Johannes H. Voigt, Bruno Bienfait, Frank Oellien, and Mare C. Nicklaus. Enhanced cactvs browser of the open nci database. Journal of Chemical Information and Computer Sciences, 42(1):46-51 57, 2002.
[81] George WA Milne, Mare C Nicklaus, John S Driscoll, Shaomeng Wang, and D Zaharevitz. National cancer institute drug information system 3d database. Journal of chemical information and computer sciences, (4(5):1219-1222 1994.
[82] Manish Monga and Edward Sausville. Developmental therapeutics program the nci: molecular target and drug discovery process. Leukemia, 16(4):520-524 2002.
[83] Amr Alhossary, Stephanus Daniel Handoko, Yuguang Mu, and Chee-Keong Kwoh. Fast, accurate, and reliable molecular docking with quickvina 2. Bioinformatics, 31(13):2214-2216 2015.
[84] Christoph Gorgulla, Andras Boeszoermenyi, Zi-Fu Wang, Patrick D Fischer, Paul W Coote, Krishma M Padmanabha Das, Yehor S Malets, Dmytro S Radchenko, Yurii S Moroz, David A Scott, et An open-source drug discovery platform enables ultra-large virtual screens, Nature, 580(7805):663-668, 2020.
[85] Christoph Gorgulla, AkshatKumar Nigam, Matt Koop, Suleyman S Cinaroglu, Christopher Secker, Mohammad Haddadnia, Abhishek Kumar, Yehor Malets, Alexander Hasson, Roni Levin-Konigsberg, et al. Virtualflow 2.0-the next generation drug discovery platform enabling adaptive screens of 69 billion molecules. bioRxiv, pages 2023-04, 2023.
[86] David Ryan Koes, Matthew P Baumgartner, and Carlos J Camacho. Lessons learned in empirical scoring with smina from the csar 2011 benchmarking exercise. Journal of chemical information and modeling, (3(8):1893-199)) 2013.
[87] Gabriel dos Passos Gomes, Robert Pollice, and Alán Aspuru-Guzik. Navigating through the maze of homogeneous catalyst design with machine learning. Trends in Chemistry, 3(2):66- 1 10, February 2021.
[88] Philippe Schwaller, Alain C. Vaucher, Ruben Laplaza, Charlotte Bunne, Andreas Krause, Clemence Corminboeuf, and Teodoro Laino. Machine intelligence for chemical reaction space. WIREs Computational Molecular Science, March 2022.
[89] Frank Jensen. Locating minima on seams of intersecting potential energy surfaces. An application to transition structure modeling. Journal of the American Chemical Society, 15665)1196-,163, February 1992.
[90] Manfred Reetz. Dyotropic rearrangements, a new class of orbital-symmetty controlled reactions. type ii. Angewandte Chemie International Edition in English, 11(2)1130-131, 1972.
[91] Amy C. Anderson. The process of structure-based drug design. Chemistry & Biology, 10(9)77777797, 2003.
[92] Xingang Peng, Shitong Luo, Jiaqi Guan, Qi Xie, Jian Peng, and Jianzhu Ma. Pocket2Mol: Efficient molecular sampling based on 3D protein pockets. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 17644-17653 PMLR, 17-23 Jul 2022.
[93] Jiaqi Guan, Xiangxin Zhou, Yuwei Yang, Yu Bao, Jian Peng, Jianzhu Ma, Qiang Liu, Liang Wang, and Quanquan Gu. DecompDiff: Diffusion models with decomposed priors for structure- based drug design. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 11827-11846. PMLR, 23-29 Jul 2023.
[94] Jiaqi Guan, Wesley Wei Qian, Xingang Peng, Yufeng Su, Jian Peng, and Jianzhu Ma. 3d equivariant diffusion for aarget-aware molecule generation and affinity prediction. In The Eleventh International Conference on Learning Representations, 2023.
16

---

[95] Vikas Nanda and Ronald L Koder. Designing artificial enzymes by intuition and computation. Nature chemistry, 2(1):15-24, 2010.
[96] Edward O. Knapp. Changwon Suh, Rafael Kömez-Bombaree Jorge Aguilera- Iparraguirre, and Alán Sapuru-Guzik, What is high-throughpu virtual screening? a perspective from organic materials discovery. Annual Review of Materials Research, 45(1):195-216) 2015.
[97] Michel Gendreau and Potvin Jean-Yves. Handbook of metaheuristics. Springer Vol. 2, 2010.
[98]
Christian Blum and Andrea Roli. Metaheuristics combinatorial optimization: Overview and conceptual comparison. ACM Comput. Surv., 35(3):268-33))) 00088-00,0000 2003.
[99] Robert Pollice, Gabriel dos Passos Gomes, Matteo Aldeghi, Riley J Hickman, Mario Krenn, Cyrille Lavigne, Michael Lindner-D' Addario, AkshatKumar Nigam, Cher Tian Ser. Zhenpeng Yao, et al. Data-driven strategies for accelerated materials design. Accounts of Chemical Research, 54(4))849-860, 2021.
[100] Daniel Schwalbe- Koda and Rafael Gómez-Bombaree Generative models for automatic chemical design. In Machine Learning Meets Quantum Physics, pages 447-467. Springer, 2020.
[101] Durk P Kingma and Prafulla Dhariwal. Glow: Generative flow with invertible 1x1 convolutions. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa- Desa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
[102] Ian Goodfellow, Jean Pouget-Abadie Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. In Advances in neural information processing systems, pages 2672-2680, 2014.
[103] Gabriel Lima Guimaraes, Benjamin Sanchez-Lengelingg Carlos Outeiral, Pedro Luis Cunha Farias, and Alán Aspuru-Guzik. Objective-reinfoceed generative adversarial networks (organ) for sequence generation models. arXiv preprint arXiv: 10843, 2017.
[104] Richard S Sutton and Andrew Barto. Reinforcement learning: An introduction. MIT press, 2018.
[105] Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural computation, 173: 1997.
[106] Alex Graves. Generating sequences with recurrent neural networks. arXiv preprint arXiv:/308.0880 2013.
[107] Junyoung Chung, Caglar Gulcehre, KyungHyun Cho, and Yoshua Bengio. Empirical evaluation gated recurrent neural networks on sequence modeling. arXiv preprint arXiv: 1412.3555, 2014.
[108] Daniel Flam-Shepherd, Kevin Zhu, and Alán Aspuru- Guzik. Language models can learn complex molecular distributions. Nature Communications, 13(1):1-10, 2022.
[109] John H Holland. Genetic algorithms. Scientific american, 267(1): 1992.
[110] Markus Hartenfeller and Gisbert Schneider. Enabling future drug discovery by de novo design.
WIREs Computational Molecular Science, 1(5):742-759, 2011.
[111] Seyedali Mirjalili. Genetic algorithm. In Evolutionary algorithms and neural networks, pages
43-55. Springer, 2019.
[112] AkshatKumar Nigam, Robert Pollice, Mario Krenn, Gabriel dos Passos Gomes, and Alan Auzuk. Beyond generative models: Superfast traversal, optimization, novelty, ex- ploration and discovery (stoned) algorithm for molecules using selfies. Chemical Science, 2021.
[113] Christoph J Brabec, Jens A Hauch, Pavel Schilinsky, and Christoph Waldauf. Production aspects of organic photovoltaics and their impact on the commercialization of devices. MRS bulletin, 30(1):0--5,, 2005.
17

---

[114] Lin X. Chen. Organic solar cells: Recent progress and challenges. ACS Energy Letters, 4(10):2537-2539, 2019.
[115] Rongming Xue, Jingwen Zhang, Yaowen Li, and Yongfang Li. Organic solar cell materials toward commercialization. Small, 14(41):180179, 2018.
[116] Carsten Deibel, Vladimir Dyakonov, and Christoph J. Brabec. Organic bulk- nullooeetrojunctionn solar cells. IEEE Journal of Selected Topics in Quantum Electronics, 16(6):11171527, 2010.
[117] Rene AJ Janssen, Jan C Hummelen, and N Serdar Sariciftei. Polymer-fullerent bulk hetero- junction solar cells. MRS bulletin, 30(1):33- 36, 2005.
[118] Brian A Gregg. The photoconversion mechanism of excitonic solar cells. MRS bulletin, 30(1):20-22, 2005.
[119] Kelly L. Famm-Ganame Richard D. Smith, James B. Dunbar, Jeanne A. Stuckey, and Heather Carlson. Csar benchmark exercise 2011-2012: Evaluation of results from docking and relative ranking of blinded congeneric series. Journal of Chemical Information and Modeling, 53(8): 1853-1870, 2013.
[120] Yan Li, Zhihai Liu, Jie Li, Li Han, Jie Liu, Zhixiong Zhao, and Renxiao Wang. Comparative assessment of scoring functions on an updated benchmark: 1. compilation of the test set. Journal of Chemical Information and Modeling, 54(6): 2014.
[121] Yan Li, Li Han, Zhihai Liu, and Renxiao Wang. Comparative assessment of scoring functions updated benchmark: 2. evaluation methods and general results. Journal of Chemical Information and Modeling, 14()):1717-1736, 2014.
[122] Zhe Wang. Huiyong Sun, Xiaojun Yao, Dan Li, Lei Xu, Youyong Li, Sheng Tian, and Tingjun Hou. Comprehensive evaluation ten docking programs on a diverse set of protein-ligand complexes: the prediction accuracy of sampling power and scoring power. Phys. Chem. Chem.
Phys., 18:12964-12975, 2016.
[123] Zoe Cournia, Bryce K. Allen, Thijs Beuming, David A. Pearlman, Brian K. Radak, and Woody Sherman. Rigorous free energy simulations in virtual screening. Journal of Chemical Information and Modeling, (009)4113-4116 2020.
[124] Tobiasz Cieplinski, Tomasz Danel, Sabina Podlewska, and Stanislaw Jastrzebski. We should at least be able to design molecules that dock well. arXiv preprint arXiv:2006. 16955, 2020.
[125] G Richard Bickerton, Gaia Paolini, Jérémy Besnard, Sorel Muresan, and Andrew L Hopkins. Quantifying the chemical beauty of drugs. Nature chemistry, 4(2):90-98, 2012.
[126] Rafael Gómez-Bombardee Jennifer N. Wei, David Duvenaud, José Miguel Hernández-Lobatt Benjamín Sänchez-Lengeling Dennis Sheberla, Jorge Hirzel, Ryan P. Adams, and Alán Aspuru-Guzik. Automatic chemical design using a data- driven continuous representation of molecules. arXiv, page 1610.02415v2, 2017.
Timothy D.
[127] Wengong Jin, Regina Barzilay, and Tommi Jaakkola. Multi-objective molecule generation using interpretable substructures. In International Conference Machine Learning, pages 4849 4859. PMLR, 2020.
[128] Christopher A Lipinski, Franco Lombardo, Beryl W Dominy, and Paul J Feeney. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Advanced drug delivery reviews, 23(1-3):3-2 25, 1997.
[129] Jerome Eberhardt, Diogo Santos-Martins, Andreas F. Tillack, and Stefano Forli. Autodock vina 1.2.0: New docking methods, expanded force field, and python bindings. Journal of Chemical Information and Modeling, 61 (8)38811-3888, 2021.
[130] Milan Voršilák, Michal Koláł, Ivan Čmelo, and Daniel Svozil. Syba: Bayesian estimation of synthetic accessibility of organic compounds. Journal of cheminformatics, 12(1):1-13, 2020.
18

---

[131] K. N. Houk and Fang Liu. Holy grails for computational organic chemistry and biochemistry. Accounts of Chemical Research, 50(3):539-544 2017.
[132] Carl Poree and Franziska Schoenebeck. A holy grail in chemistry: Computational catalyst design: Feasible or fiction? Accounts of Chemical Research, 50(3):605-608, 2017.
[133] Sharon Hammes- Catalysts by design: The power theory. Accounts of Chemical Research, (063):561-5500 2017.
[134] Todd Martínez. Ab initio reactive computer aided molecular design. Accounts of Chemical Research, 50(3):652-66) 656, 2017.
[135] Keith J. Laidler and M. Christine King. Development of transition-state theory. The Journal of Physical Chemistry, 87(15):2657-2664 1983.
[136] Donald G. Truhlar, Bruce C. Garrett, and Stephen J. Klippenstein. Current status of transition- state theory. The Journal of Physical Chemistry, (0031))12771-122800 1996.
[137] A. D. McNaught, A. Wilkinson, and S. J. Chalk. IUPAC. Compendium of Chemical Terminol- ogy, (the "Gold Book"), Online Version. Blackwell Scientific Publications, Oxford, second edition, 2019.
[138] Pascal Friederich, Gabriel dos Passos Gomes, Riccardo De Bin, Alán Aspuru-Guzio and David Balcells. Machine learning dihydrogen activation in the chemical space surrounding Vaska's complex. Chemical Science, 1111))4554 4601, 2020.
[139] Colin A. Grambow, Adeel Jamal, Yi-Pei Li, William H. Green, Judit Zádor, and Yury V. Suleimanov. Unimolecular reaction pathways of a y-ketohydroperoxide from combined application of automated reaction discovery methods. Journal of the American Chemical Society, 14053):0035-1048 2018.
[140] Lagnajit Pattanaik, John B. Ingraham, Colin A. Grambow, and William H. Green, Generating transition states of isomerization reactions with deep learning. Physical Chemistry Chemical Physics, 22(41):23618-22 226226, 2020.
[141] Kjell Jorner, Tore Brinck, Per-Ola Norrby, and David Buttar. Machine learning meets mech- anistic modelling for accurate prediction of experimental activation energies. Chem. Sci., 12:1163-1175, 2021.
[142] Yunhan Chu, Wouter Heyndrickx, Giovanni Occhipinti, Vidar R. Jensen, and Bjøm K. Alsberg. An evolutionary algorithm for de novo optimization of functional transition metal compounds. Journal of the American Chemical Sociery, 134(21):8885 8895,, May 2012.
[143] Ruben Laplaza, Simone Gallarati, and Clemence Corminboeuf. Genetic optimization of homogeneous catalysts. Chenhisryy-Mehhods, 2(6), June 2022.
[144] Julius Seumer and Jan Jensen. Computational evolution of new catalysts for the morita- billinn reaction. ChemRxiv, 2022.
[145] J. N. Bronsted and K. Pedersen. Die katalytische zersetzung des nitramids und ihre physikalisch-chemische bedeutung. Zeitschrift für Phys. Chemie, Stöchiometrie und Ver- wandtschaftslehre, 108:185-235, 1924.
[146] Ronald Percy Bell. The theory of reactions involving proton transfers. Proceedings of the Royal Society of London. Series Mathhematccal and Physical Sciences, 154(882)4444429, 1936.
[147] M. G. Evans and M. Polanyi. Further considerations the thermodynamics of chemical equilibria and reaction rates. Trans. Faraday Soc., 32:1333-1360 1936.
[148] Joseph R. Murdoch. A simple relationship between empirical theories for predicting barrier heights of electron-, proton-, atom-, and group-transte reactions. Journal of the American Chemical Society, 105(8):155-2266 April 1983.
19

---

[149] Wengong Jin, Regina Barzilay, and Tommi Jaakkola. Hierarchical generation of molecular graphs using structural motifs. In International conference on machine learning, pages PMLR, 2020.
[150] David Rogers and Mathew Hahn. extended-connecctvity fingerprints. Journal of chemical information and modeling, 50(5):742-754, 2010.
[151] AkshatKumar Nigam, Robert Pollice, Matthew F. D. Hurley, Riley Hickman, Matteo Aldeghi, Naruki Yoshikawa, Seyone Chithrananda, Vincent A. Voelz, and Alán Aspuru-Guzik, Assigning confidence to molecular property prediction. Expert Opinion on Drug Discovery, 1609):100--1023, 2021.
[152] Thomas J. Penfold. On predicting the excited-state properties of thermally activated delayed fluorescence emitters. The Journal of Physical Chemistry C. 11922)):353551333 2015.
[153] Youichi Tsuchiya, Keita Tsuji, Ko Inada, Fatima Bencheikh, Yan Geng, H. Shaun Kwak, Thomas J. L. Mustard, Mathew D. Halls, Hajime Nakanotani, and Chihaya Adachi. Molecular design based on donor-weak donor scaffold for blue thermally-activated delayed fluorescence designed by combinatorial dft calculations. Frontiers in Chemistry, 8, 2020.
[154] Israel Fernández, Fernando P. Cossío, and Miguel A. Sierra. Dyotropic reactions: Mechanisms and synthetic applications. Chemical Reviews, (0((1))6687-6711 2009. PMID: 19746971.
[155] Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In ICLR (Poster), 2015.
[156] Markus C Scharber, David Mühlbacher, Markus Koppe, Patrick Denk, Christoph Waldauf, Alan Heeger, and Christoph Brabec. Design rules for donors in bulk-heterojunction solar towards 10% energy-conversion efficiency. Advanced materials, ((6))789-794) 2006.
[157] Noel M O' O'Boyle, Michael Banck, Craig A James, Chris Morley, Tim Vandermeersch, and Geoffrey R Hutchison. Open Babel: An open chemical toolbox. Journal of Cheminformatics, 3(1):33, December 2011.
[158] Philipp Pracht, Fabian Bohle, and Stefan Grimme. Automated exploration of the low-energy chemical space with fast quantum chemical methods. Physical Chemistry Chemical Physics, 22(14)7119-7199, 2020.
[159] Sebastian Spicher and Stefan Grimme. Robust atomistic modeling of materials, organometallic, and biochemical systems. Angewandte Chemie International Edition, 59(36):15655115755 September 2020.
[160] Henri Theil. A rank-invarinn invariant method of linear and polynomial regression analysis. Indaga- tiones mathematicae, 12(85):173, 1950.
[161] Pranab Kumar Sen. Estimates of the regression coefficient based on kendall's 's tau. Journal of the American Statistical Association, 63(324):1379 1968.
[162] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg. J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825-288 28300, 2011.
[163] Greg Landrum et al. Rdkit: Open-source cheminformatics, 2006.
[164] Qiming Sun. Libeint: An efficient general integral library for gaussian basis functions. Journal of Computational Chemistry, 36(22): 164-16611 2015.
[165] Qiming Sun, Timothy C. Berkelbach, Nick Blunt, George H. Booth, Sheng Guo, Zhendong Li, Junzi Liu, James D. McClain, Elvira R. Sayfutyarova, Sandeep Sharma, Sebastian Wouters, and Garnet Kin-Lie Chan. Pysef: the pashon-based simulations of chemistry framework. WIREs Computational Molecular Science, 8 3(1)ee3340 2018.
20

---

[166]
Qiming Sun, Xing Zhang, Samragni Banerjee, Peng Bao, Mare Barbry, Nick S. Blunt, Niko- lay A. Bogdanov, George H. Booth, Jia Chen, Guo, Jan Hermann, Matthew R. Hermes, Kevin Koh, Peter Koval, Susi Lehtola, Zhendong Li, Junzi Liu, Narbe Mardirossian, James D. McClain, Mario Motta, Bastien Mussard, Hung Q- Pham, Artem Pulkin, Wirawan Purwanto, Paul J. Robinson, Enrico Ronca, Elvira R. Sayfut- yarova, Maximilian Scheurer, Henry F. Schurkus, James E. Smith, Chong Sun, Shi-Ning Sun, Shiv Upadhyay, Lucas K. Wagner, Xiao Wang, Alec White, James Daniel Whitfield, Mark J. Williamson, Sebastian Wouters, Jun Yang. Jason M. Yu, Tianyu Zhu, Timothy C. Berkelbach, Sandeep Sharma, Alexander Yu. Sokolov, and Gamet Kin-Lic Chan. Recent developments in the pysef program package. The Journal of Chemical Physics, 153(2)0024109 2020.
Cui, Janus Eriksen, Yang Gao, Sheng
[167] P. Jeffrey Hay and Willard R. Wadt. Ab initio effective core potentials for molecular calcula- tions. potentials for k to au including the outermost core orbitals. The Journal of Chemical Physics, 82(1):299-310, 1985.
[168] G Madhavi Sastry, Matvey Adzhigirey, Tyler Day, Ramakrishna Annabhimoju, and Woody Sherman. Protein and ligand preparation: parameters, protocols, and influence on virtual screening enrichments. Journal of computer-aice molecular design, 27(3):221-255 2344, 2013.
[169] Noel M O'Boyle, Michael Banck, Craig A James, Chris Morley, Tim Vandermeersch, and Geoffrey R Hutchison. Open babel: An open chemical toolbox. Journal of cheminformatics, 3(1):1-14, 2011.
[170] Sunghwan Kim, Paul A Thiessen, Evan Bolton, Jie Chen, Gang Fu, Asta Gindulyte, Lianyi Han, Jane He, Siqian He, Benjamin A Shoemaker, al. Pubchem substance and compound databases. Nucleic acids research, 44(D1):11202-D12211 2016.
[171] Thomas A. Halgren. Merck molecular force field. I. Basis, form, scope, paramcterization, and performance of MMFF94. Journal of Computational Chemistry, 17(5-6):400-51, April 1996.
[172] Xiaolei Zhu, Keiran C. Thompson, and Todd J. Martínez. Geodesic interpolation for reaction pathways. The Journal of Chemical Physics, 150(16):16413,, April 2019.
[173] Lee-Ping Wang and Chenchen Song. Geometry optimization made simple with translation and rotation coordinates. The Journal of Chemical Physics, 144(21)214108, June 2016.
[174] F. Pedregosa, G. Varoquaux, A. Gramfort, Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg. J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825-2830, 2011.
21

---

Supplementary Information:
TARTARUS: Practical and Realistic Benchmarks for Inverse Molecular Design
S1
Introduction
Traditionally, property-guided optimization has relied on expert intuition 995 and several rounds of trial, error, and human- -nnsiird optimization, occasionally supported by computer simulations. Alternatively, computer- aasputer-assisted approaches have employed virtual screening [96] optimization algorithms such as genetic algorithms (GAs) 197, 98 [5]. More recently, with the surge of deep learning, deep generative models have emerged, specifically designed to operate in chemical space and tackle inverse molecular design [99 100) 1. This has led to the development of numerous algorithmic approaches for this purpose, with the most popular including variational autoencoders (VAEs) [1] generative adversarial networks (GANs) [02 103], and reinforcement learning (RL) 104 12
S2
Methods Overview
In this section, we provide an overview of the molecular generative models employed throughout this work and summarize the associated design choices needed to make during their implementation. The molccular design algorithms we considered are VAEs, long short-term memory hill climbing
(LSTM-HC) models [O05] 33 31, REINVENT 30, JANUS 77. and graph-based genetic algorithm
(GB-GA) 34 At the core of the majority of these approaches are molecular string representations, the most commonly used of which is the Simplified Molecular Input Line Entry System (SMILES) [46]. Accordingly, many of the algorithms tested rely on predicting subsequent characters from partial strings to propose structures. However, algorithms based on SMILES can regularly produce invalid strings that do not represent molecules, which is problematic both in terms of robustness and interpretability of the corresponding methodologies 19 20 Consequently, this issue was addressed systematically by introducing Self- Self-Referencing Embedded Strings (SELFIES) [19], molecular string representation that guarantees validity. Thus, unlike for SMILES, every arbitrary combination of SELFIES characters represents a molecule. Nevertheless, its impact on structure optimization has not yet been evaluated systematically (0)). To address this issue, we modify some of the existing generative models relying on SMILES to be also compatible with SELFIES and test their performance depending on representation, similar to how has been done recently [8]
Among the models tested, REINVENT, the VAEs, and the LSTM-HC models use recurrent neural networks (RNNs) [06] to learn the conditional probability distributions of the characters that represent molecules. RNNs are a class of artificial neural networks (ANNs) that utilize sequential information from their previous predictions and states. They have been incorporated into molecular design algorithms with special NN node architectures such as gated recurrent units (GRUs) [077] and LSTM cells (105]. The first of these models, i.c., REINVENT [30), is an RL-based approach that relies on an LSTM- RNN as an agent which is tasked with generating molecules with desired properties. Over time and continued training, the agent learns to propose compounds with increasingly favorable target property values.
For the VAEs, used the implementation described by K6meez-Bobarrell et al. starting point (3I) Therein, SMILES are converted to one-hot encodings and, subsequently, passed through I-dimensional convolutional neural network (CNN) encoder that generates a continuous latent space. A separate RNN decoder with GRU nodes regenerates the SMILES strings from the latent space. Molecular optimization is performed on the continuous representation of the latent space in stochastic and iterative manner. Specifically, the best available molecule for given task is encoded into the latent space of trained VAE. Subsequently, random Gaussian noise is added to the obtained latent vector to produce population of latent vectors that can be decoded to a population of new candidate structures. These structures are evaluated and the best compound obtained used as a seed in the subsequent iteration. In our experience, this strategy leads to more stable optimization compared to direct gradient-based optimization in the latent space as described in the literature B1]. Importantly, we modified the original implementation to be compatible with both SMILES
S1

---

and SELFIES (cf. Computational Details in the Supporting Information). These variations will be referred to as SMILES- VAE and SELFIES-VAL respectively.
rely on LSTM cells in the NN architecture to model the character sequence probability distributions. After initial training, the resulting model is sampled using randomly truncated strings of the best currently available molecules for given task. These truncated strings are used as seeds and completed stochastically by sampling the learned conditional probability distributions to produce a population of candidate structures. The best new compounds obtained after property evaluation are used to retrain the model and initiate a subsequent sampling cycle. Thus, iterative sampling and retraining gradually improves the designs of the LSTM-HC models. Again, we test both the original implementation of LSTM-HC that relies on and modified version making Notably, the AMILES-LSTMAA model was
As their name implies, the LSTM-HC models 105
3
T08
SMILES 33 3 which will be referred to as use of SELFIES, referred to as SELFIES-LSTM-HHE among the best performing models in the original publication of the GuacaMol benchmarks (referred to as "SMILES LSTM" in that reference).
Before the increasing adoption of ML approaches for molecular design, GAs [009] among the most popular computer- based methods for that purpose. Inspired by natural selection as defined in Darwinian evolution, GAs are heuristic population-base rely on repeated stochastic generation of offspring from an existing population of candidate solutions via the genetic operations mutation and crossover. Subsequently, selection of the candidate solutions with the highest fitness determines the population to be propagated to the next generation, starting another iteration. In this work, we consider two specific implementations for inverse molecular design, GB-GA [34), and JANUS 77]. The former leverages structural information from a dataset of reference molecules by mimicking the corresponding distribution of atoms and bonds when performing genetic operations. Notably, in the original publication of the GuacaMol benchmarks [3) GB-GA achieved the best overall performance. In contrast, JANUS is GA augmented by an NN classifier (7] (referred to as "JANUS+C" in the original publication) that actively judges the quality of molecule before performing full fitness evaluation. Molecules classified as likely possessing high fitness are then propagated to the next generation for subsequent evaluation. For mutation and crossover, JANUS relies STONED-SELFIEE 1112), of algorithms based on the guaranteed validity of SELFIES to perform structural modifications. Additionally, unlike GB-GA, JANUS maintains two distinct molecule populations, one explorative and exploitative with respect to conducted structural changes 77. These populations exchange some members at every iteration but otherwise explore the chemical space independently
were
optimization algorithms. They
When using TARTARUS, the following procedures should be adopted to obtain benchmark results that are consistent with the ones provided herein. The first step for running one of the benchmarks, necessary, is to train the generative model on the provided dataset. For all the ML models, we used the first 80% of the reference molecules for training and the remaining 20% for hyperparameter optimization. Then, the (trained) model is tasked with proposing structures to be evaluated by the objective function of the corresponding benchmark task. Notably, structure optimization was always initiated using the best reference molecule from the corresponding dataset. For the benchmarks concerned with designing photovoltaics, organic emitters, and protein ligands, structure optimization was carried out with a population size of 500 and limit of 10 iterations, leading maximum number of 5,000 proposed compounds overall. For the design of chemical reaction substrates, we used the same maximum number of proposed compounds but used a population size of 100 and limited the number of iterations to 50 instead. Additionally. the associated run time was limited to 24 hours, which resulted in termination for several molecular design runs before reaching 5,000 molecule evaluations. Furthermore, to increase robustness and reproducibility of our results, we repeated each optimization run five times, allowing us to report the corresponding outcomes with both an average and a standard deviation. We believe that this pesource-constraie is necessary for fairly comparing methods and should be used as a standard by the community, A detailed account of the parameters and settings used for running each of the models is provided in the Computational Details section of the Supporting Information.
comparison approach
S2

---

S3
Additional Results
S3.1 Design of Organic Photovoltaics
OPVs offer simplified, cost-effective production 113 35 and enhanced mechanical properties,
particularly regarding specific weight and flexibility still have lower power conversion efficiencies (PCEs) and shorter device lifetimes [144) prompting
38
39]. Despite significant progress, OPVs
ongoing research efforts in molecular design for OSCs 115
114 39
The simplest OPV cells comprise two electrodes and photoactive layer for photoconversion, typically containing two distinct materials: donor and acceptor [35]. bulk heterojunction cells [16]. donors and acceptors are mixed, allowing nanoscale phase separation to maximize interfacial area and minimize distance within the phases [17]. Upon light exposure, excitons-neuttee in the photoactive layer with limited lifetimes and diffusion lengths [117]. The bulk heterojunction architecture enables excitons to reach the interface, dissociating into a freee electron and hole that become charge carriers across the phases [18 117 35 These charge carriers travel through the photoactive layer to the electrodes, generating photocurrent [55. For exciton charge separation the donor-acceptor interface to occur, the process must be energetically neutral or favorable [188 (35]. Therefore, the exciton energy, estimated by the HOMO-LUMO gap of the light-absorbing material, must be greater than or equal to the effective heterojunction bandgap [118 35 approximated by the energy difference between the donor's highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital of the acceptor (LUMO) 41. Notably, PCE is the percentage of power from the incident solar irradiation that is converted into electricity by the OPV device. We believe that the reduced reference dataset presents a more realistic scenario for new molecular design projects and is better suited for benchmarking generative models, especially when property simulations are time-consuming. We propose the following two molecular design benchmark objectives:
quasi-particles consisting of bound electron-hole pairs [188 351 form
1. Maximize the following function: PCEPсBM SAscore.
2. Maximize the following function: PCFPCdтBT Ascore.
GB-GA
JANUS
N
N
N
PCEpcaM" SAscore 8.056
PCE SAscore 7.851
SMILES-LSTM-HC
JANUS
PCEpcoTer SAscore 32.057
PCE SAscore 31.717
Figure S1: Best molecules found each of the benchmark tasks inspired by the design of organic photovoltaics. Additionally, the corresponding objective values and the molecular design models that proposed the structures are indicated.
S3

---

The best molecules found in each of the design of organic photovoltaics benchmark tasks together with the corresponding objective values and the model that proposed them are depicted in Figure ST
S3.2 Design of Organic Emitters
The next set of benchmarks is inspired by the design of purely organic emissive materials for organic light-emitting diodes (OLEDs), which received significant attention in recent years [55 after the discovery of thermally activated delayed fluorescence (TADF) in the field [60]. Their main applications digital screens and lighting devices [55]. For the former application, compared to alternative technologies. OLEDs offer improved image quality and enable both lighter and thinner devices [55]. For the latter, OLEDs are potentially more energy-efficient [55]. In OLED devices relying on TADF, after electric excitation, light generation takes place from the first excited singlet state via fluorescence (60). However, only 25 % of the initial excitons are excited via an electric current into singlet states and contribute directly to light emission via excited state thermal relaxation and subsequent prompt fluorescence (55 60) The 75 % of the initial excitons that excited into triplet states relax thermally to the corresponding first excited triplet states [55 60). However, in ordinary organic molecules, light emission from their first excited triplet states via phosphorescence is slow, giving rise to various radiationless decay and decomposition pathways 57 Consequently, both device efficiency and device lifetime are reduced considerably unless these triplet excitons can still be utilized for light production [55) To achieve that, TADF emitters rely on minimizing the energy difference between the first excited singlet and triplet states, i.c., the singlet-triplet gap, allowing thermal upconversion of the triplet excitons to the first excited singlet state giving rise to delayed fluorescence 57 Importantly, under the assumption of fast both forward and reverse intersystem crossing (ISC), the steady-state triplet population is governed by the singlet-triplet gap and its reduction leads to reduced triplet population and acceleration of delayed fluorescence. This increases internal quantum efficiency up to a maximum of 100 % and reduces decomposition of the emissive material (55)57]. Designing efficient emissive organic materials for blue OLEDs is of particular interest. This can be achieved by targeting organic molecules with excitation energies between their ground state and their first excited singlet states that correspond to tasks are summarized as follows:
the energy of blue light. The fitness functions of these three
Minimize the singlet- gap, AE(ST-Ti).
Maximize the oscillator strength for the transition between S and So. f12-
Maximize the following function:
+f12 - AS- - ( S1) - 3.2 eV
The best molecules found in each of the design of organic emitters benchmark tasks together with the corresponding objective values and the model that proposed them are depicted in Figure S2
S3.3 Design of Protein Ligands
Notably, while docking simulations are standard tool in virtual screening pipelines of drug discovery campaigns, their accuracy compared to experimental binding affinities is at best modest 122. Thus, typical workflows use them merely for a preselection which narrowed down further with subsequent free energy simulations [13]]. Nevertheless, for molecular design benchmarks, docking still provides the best trade-off between computational efficiency and relevance to real-world molecular design. Additionally, it is important to realize that drug design requires the consideration of many other molecular design aspects that make or break a hit compound, e.g.. toxicity, solubility, stability and many more. However, these properties model computationally, we decided to disregard them for our set of benchmarks. Notably, finding small molecule ligands for the selected proteins marks the first step towards the development of treatments for various important conditions.
significantly harder to
Importantly, most likely due to its maturity, complexity, and demand for resources, drug design was probably the first chemical problem where molecular design algorithms were tested comprehensively.
S4

---

JANUS
NH2
NH
NH
Singles triplet 0. 000
Singiet trppet gap 0,012 ev
JANUS
CH2
Osciliator Streegth - 2.45
Csciliaoo Strongth = 2.45
OH
HD
HN
Combined Objective 0.18
Combinad Objective 0.07
Figure S2: Best molecules found each of the benchmark tasks inspired by the design organic emitters. Additionally, the corresponding objective values and the molecular design models that proposed the structures are indicated.
The use of computer algorithms has a long-standing history in medicinal chemistry and GA-based molecular design algorithms making use of full atomic representations have already been used as early as 1993 [6]. Initial toy tasks for testing these algorithms included rediscovery of known drugs via structural similarity metric, highly simplified molecular docking of ligands to rigid protein binding sites minimizing the interaction scores, and the optimization of estimated molecular properties like the decadic logarithm of the t-octanol-watee distribution coefficient (log P) 66. Interestingly, some of the still most widely used benchmarks for generative models rely on essentially the same types of metries as they are simple to implement and compute 1,
3 124 29). More recently, the quantitative estimate of drug-likeness (QED) was introduced, which is a desirability function that uses the position of a small set of common and simple molecular descriptors relative the corresponding distributions in dataset of approved drugs to estimate structural resemblance to therapeutics [125]. Due to its simplicity, it found immediate application in several molecular design benchmarks [126) 103 [27]. However, using QED alone in generative models is not meaningful for finding drug candidates as it only accounts for the general structural requirements of drug-like molecules but disregards intended modes of action with respect to specific targets entirely. The specific benchmark objectives we implemented are summarized below.
Minimize the docking score to ISYH, ΔEIYY
Minimize the docking score to 6Y2F, ΔEGYFF.
Minimize the docking score to 4LDE, ΛEADDE
S5

---

Notably, the corresponding objective functions do not solely consist of the docking scores but also have hard structural constraints directly incorporated. When these constraints are not fulfilled, an extremely unfavorable score of 10,000 is returned instead of the actual docking score. They consist of a set of filters that checks for the presence of unstable or reactive structural moieties and determines whether the compound in question fulfils Lipinski's Rule of Five [228]. Notably, most of these filters were developed based on our previous experience using molecular design algorithms to minimize docking scores and are tailored to avoid extremely unstable and reactive molecules that seem to be strongly favored by molecular docking simulations 77]. Additionally, the constraints avoid rings with more than 8 members the docking approach implemented is unable to sample the corresponding conformations proper manner [19]] To fulfill them, the proposed structure needs to have an SAscore smaller than 4.5, which is the revised optimal threshold for that metric proposed in the literature [30), and a QED value larger than 0.3, which corresponds to the first quartile of the distribution of QED values for compounds in the CHEMBL database 125 A list of these metrics is provided here:
Absence of reactive groups.
Absence of formal charges.
Absence of radicals.
At most 2 bridgehead atoms.
No rings larger than 8-membered.
Fulfils Lipinski's Rule of Five.
S Ascore < 4.5.
QED > 0.3.
TPSA> 140.
Molecule passes the PAINS and WEHI and MCF filters.
Molecule does not contain Si and Sn atoms.
The best molecules found in each of the design of protein ligands benchmark tasks together with the corresponding objective values and the model that proposed them are depicted in Figure S3
$3.4 Design of Chemical Reaction Substrates
Whereas, classically, the optimization of reaction parameters was largely dominated by experimental work, in recent years, the significant increase in computing power and the continuous improvement of computer algorithms enabled molccular simulations to play an increasingly important part [33)
With the aid of transition state (TS) theory, fundamental reaction parameters such as thermodynamic feasibility, reaction rate and selectivity can be computed from first principles [35]] 136). This requires explicit modeling of the corresponding TS, postulated state the multi-dimensionll potential energy surface (PES) of the process, which lies on a saddle point of order one [37]. Due to the
difficulty of finding such saddle points in high-dimensional spaces and the often delicate electronic structure associated with the corresponding structures, in practice, automated TS well-behaved case studies combined with robust methodologies still have some room for improvement in that
often suffer from high failure rates in the range of 10-50% 140) and optimizations
respect[[4]]. Additionally, they are typically carried out with relatively resource-intensive density functional approximation (DFA) calculations taking on the order of hours or even days to complete [38 140.. Overall, these issues make them ill-suited for benchmarking molecular design algorithms or for routine application combined with generative models in computer- cenputer-guided guided inverse molecular design campaigns. Nevertheless, GAs have been employed for computational catalyst design, particularly via the use of regression models based electronic structure descriptors derived from SQC simulations 142 and based on both structural and electronic structure descriptors obtained from DFA computations [43]. Most notably, very recently, GAs have been employed to optimize an organocatalyst for the Morita-Baylis-Ii Hilman reaction 144
S6

---

1SYH
Qutiiiine 2 Seare: keaimsl
Suore
Smina Score
5-12.0 koaimol
JANUS CH2
HJO
6Y2F
Hot
HLC
"IrcH
HjN
Qankiiis Score 11 lecattre Smera Score
CuickV/r 2 Searw -11 5 koaiiel Biina
:-12 kealtroll
GB-GA
RERIVENT
CH+
H2o
4LDE
DH
CuéllV
Figure S3: Best molecules found each of the benchmark tasks inspired by the design of protein ligands. Additionally, the corresponding objective values and the molecular design models that proposed the structures are indicated.
The two force fields cross at approximately the TS geometry, and we introduce a coupling term that turns this into avoided crossing, ensuring smoothness of the PES [89]. The resulting PES has local maximum on the ground state surface (the TS) and corresponding local minimum on the excited state surface, at approximately the same geometry. This allows optimizing TS geometries with robust gradient-basee activation energy estimates for reasonably sized molecules within minutes. and, in our hands, reaches very high success rates above 99.9% on set of test reactions. Notably, we implemented this method our Python package polanyi, which will be described in more detail in separate publication. To generate the reference dataset for this set of benchmarks, starting from the unsubstituted reactive core structure, we performed repetitive cycles of STONED-SEE -ELFIIS mutations [112] followed by removing all proposed compounds that violated the core and functional group constraints (Details the Supporting Information). Thus, after several cycles, we obtained approximately 60,000 molecules defining the reference structures, which we refer to SNB-60K dataset. Notably, in the selected reaction, there is only one TS connecting reactants and products in the selected reaction. Additionally, the fourth task aims to break the Bell-Evans -Polnyi principle [145)148), linear free energy relationship that holds empirically for large number of reactions. The following list summarizes the four benchmark tasks for chemical reactivity.
minimization algorithms. Thus, the SEAM force field method delivers
Minimize the activation energy, AE and maintain the and functional group constraints.
Minimize the reaction energy, A E,,, and maintain the core and functional group constraints.
Minimize the following function:
+AE + AE,, and maintain the core, functional group and SAscore constraints.
Minimize the following function:
-AE + AE, and maintain the core, functional group and SAscore constraints.
S7

---

On top of these primary objectives, we added several hard constraints that the target molecules need to fulfiil in order to reward generative models that propose realistic and feasible molecules. In particular, all substrates need to retain the syw-sesquinorborrene motif (referred to as "core constraint") which required for the reaction to take place. Additionally, we selected a set of unstable and reactive substructures that need to be avoided (referred to "functional group constraint", details in section $5.2.4 of the Supporting Information). Furthermore, for the two benchmark tasks with objective functions combining two target properties, we also required all proposed structures to possess an SAscore 4d of 6.0 or lower (referred to as "SAscore constraint"). The following list summarizes the four benchmark tasks for chemical reactivity.
JANUS
GB-GA
OH
NHz
AE' 44.04 kca//mol
AE' 51.28 kcalmol
JANUS
GB-GA
HO
SH
NH2
HN
NH2
AE, = -9. 10 kcal/mol
AE, -48 76 kcalmmol
JANUS
GB-GA
NH
HO
+4E' 32.79 kcal/mol
+4E' + 35. 47 kcal/mol
GB-GA
JANUS
HN
NH2
-AE' + AE, -10..51 koal/mol
-AE' AE, = 98.81 kcal/mol
Figure S4: Best molecules found in each of the benchmark tasks inspired by the design of chemical reaction substrates. Additionally, the corresponding objective values and the molecular design models that proposed the structures are indicated.
The best molecules found in each of the design of chemical reaction substrates benchmark tasks together with the corresponding objective values and the model that proposed them are depicted in Figure S3
S8

---

$3.5 Model Timing Comparison
For ML-based molecular design algorithms, the pre- conditioning corresponds to training time and sampling time, respectively. For the GA-based algorithms considered, they translate to the duration of the derivation of genetic operators from a reference dataset, and of applying the genetic operators to propose new candidate solutions, respectively. When the number of property evaluations is kept constant, assuming that the molecular size distribution between the generative models does not differ significantly, these steps the major origin of timing differences between the algorithms. Notably, as was done for generating the molecular design results, we derived these timing metrics from five independent measurements and provide the results averages with standard deviations.
Comparison of the single epoch times with and without GPU allows estimating which models profit most from GPU usage. find that the longer the single epoch time, the more the model profits from using a GPU which is consistent with expectations. Finally, looking at the sampling times, we find that REINVENT significantly outperforms all other methods considered needing less than minute. Most of the other molecular design algorithms need between 2 and 3 minutes with GB-GA having the slowest sampling time of 6 minutes. Nevertheless, the differences between the methods are less pronounced here.
Table S5: Raw values for the timing benchmarks. Mean and standard deviation (mean + s.d) timings for different models are provided based on five independent runs. N.A. means not applicable.
TRAINING TIME (S]
EPOCHS 75.6
SAMPLE TIME (s)
CPU EPOCH TIME [8]
GPU EPOCH TIME (8)
SEE FIRS-VAV SMILES- FA-
S6910 912
5.9
301
5
10810 949
538 36
34868 667
24.1 1..1+t00 6.2
154
25476 724
$15 + 59
MOFLOW
2804 + 216
70.1 5.4
62.41 =
2131 63
45 +
GB-GA
N.A.
N.A.
314
3.653 0.007
N.A.
JANUS
N.A.
N.A.
147 2:4
10.3 2.4
N.A.
RRINVENT
2844
1.2
81.7
STM-RC
10.5
139
FIES-LSTM-HC
Estimated CPU Training Time
Time (seconds)
30a
30"
Figure S5: Estimated CPU training time the deep generative models based on the CPU single epoch time and the number of epochs during training with GPUs. Models are trained on subset of DTP Open Compound Collection. Results are provided as mean (main bar) from five independent The numbers the abscissa each refer to the following molecular design algorithms. Model SELFIES- VAE; Model 1: SMILES-VEL VAE; Model 2: REINVENT; Model 3: SMILES-LSTM-HCC Model 4: SELFIES-LSTM-HC.
The numerical results of all the timing benchmarks are provided in Table S5 estimation of total CPU training times illustrated in Figure $5
The results of the
S3.6 Diversity Calculation
Following the definition of diversity in the literature 1449, for each task, we calculate diversity the proposed moleucles using the following equation:
S9

---

2 Diversity = nnn - 1) X.Y
sim X, Y)
The expression sim(X, Y) computes the pairwise molecular similarity for all structures calculated as the Tanimoto distance of the Morgan fingerprints, which were obtained with a radius of size 3 and 2048 bit size [150]. The individual results for the diversity evaluations of all the benchmark tasks are provided in Tables S6 59 Overall, we observe mostly relatively small differences in the diversity of the proposed molecules between the different models, except for the chemical reaction substrate design benchmarks. We hypothesize that the strict structural requirement of the corresponding tasks responsible for this observation only molecules with the correct reactive core structure subjected to the property simulation workflow. Additionally, we observe that all generative models except for SMILES-LSTME
show consistently a relatively high diversity across all four groups of benchmarks. The molecules proposed by AMILES-LSTME for the chemical reaction substrate benchmarks have quite low diversity compared the other generative models. Nevertheless, these results suggest that the diversity metric can be insightful in cases of especially low values, but otherwise does not allow to distinguish between the performance of the generative models in a reliable way.
for the organic photovoltaics and
Table S6: Diversity results for the organic photovoltaics benchmarks. Models are trained on a subset of the Harvard Clean Energy Project Database. Results are provided as mean and standard deviation of the best target objective values that are obtained in five independent runs in the form mean = standard deviation. Metries: РCЕP' PCBM power conversion efficiency; PCEPсотвт: PCDTBT power conversion efficiency; SAscore: synthetic accessibility score.
Pl K.bore
PC Eecurr Selnore
MARLOW
86.9 5.2%
55.0 2:24
SMI
11-51:7333 7.3%
REINVENT
OB-GA
85.12004
89.1 0.5%
Table S7: Diversity results for the organic emitters design benchmark objectives. Models are trained on subset of the GDB-13 dataset. Results are provided as mean and standard deviation of the best objective values from five independent runs in the form mean + standard deviation. Metrics: AE Ti): singlet-triplet gap; fi2: S1 and So transition oscillator strength; AE(So-S)): So and Si vertical excitation energy.
1/n AED SAISS
SELIEES VAL
950 0.2%
$9.5 1 1.6%
92.1 0.85
MORLOW
93.2 1.4%
39.1 1 0...5.0....
RE. 0.0.
SMILDS-LSTM-IC
80.91
REINVENT
00.0000.5
1.111
GB GA JANUS
01.020033
81:1::0:00 1:1:::0:0
03.1 0.45
9..1100.1
93.6 0.15
0.29
Table S8: Diversity metrics for protein-ligand design benchmarks, based models trained on subset the DTP Open Compound Collection. Metrics show mean and standard deviation of optimal target objective values over five independent runs (mean + standard deviation). Benchmark metrics: AEx denotes docking score for protein target X.
SErsy
Aovsr
SMILES VAB SELFIES VAE
0.2%
1..4%
92.8 0.4%
91.9 4: 0.7%
90.6
MoFLOW
92.5
92.5 1.5%
92.6 1.6%
SMILES -ETLE-C-LTM-CC
91.191.0.013
91.14 1 0.04%
91.15 0.002
SELFIES-LSTM-RC
01.23330.00015
REINVENT
9.22
92.4
0
G8-GA JANUS
0.3%
92.2 + 0.2%
92.2 0.1%
91.5 0.3%
91.5 + 0.7%
90.9 + 1.0%
S10

---

Table S9: Diversity results for chemical reaction substrate design benchmarks. Models, trained on benchmark dataset generated with BONED--ELFFIE, cycles, yield mean and standard deviation of optimal objective values over five runs (mean + standard deviation). Benchmark metrics: AE Activation energy of the reaction; A E,: Reaction energy.
AE
AE,
AE + AE.
-A AE,
SMILES-VAE
70.4 10.0%
81.7 1.6%
66.3 10.321 1.3%
64. = 5.13
SELFIES-VAS VAE
69.2
5.9% 66.9
64.3 4.3%
68.4
12.6%
MOFLOW MMILES-LSTM-HG
90.6
84,9 4.4%
85.7 = 3.7%
34.2
37.6 11.2%
3.8 1.5%
18.4 = 15.4%
SELFIES LSTME---LTM-II HC
59.3 11.3%
55.6 19.8%
72.5 1 5.8%
REINVENT
88.0 0.4%
88.3 0.3%
83.2 0.8%
83.2 0.8%
GB-GA
1.1%
83.3 2.2%
76.3 1.6%
80.8 + 0.6%
JANUS
77.3 4: 2.6%
81.2
1.4%
68.9 4: 4.0%
71.4 + 2.0%
S4 Additional Discussion
In this section, we will put our findings into perspective and compare them to the pre-existing knowledge in the field. Specifically, will address two major talking points. First, we will compare the molecular design benchmarks developed as part of TARTARUS to existing benchmarks that are currently being used. Second, we will investigate the results obtained in the various molccular design objectives in detail and discuss their potential implications for the current status of artificial molecular design and future directions for the field.
S4.1 Benchmarking with TARTARUS
First, we would like to emphasize that the benchmark results provided as part of this work largely serve as a demonstration of how the objectives can be used to compare the performance of different algorithmic approaches. They should not be viewed as final performance judgement for any of the methods used as model hyperparameters were not optimized comprchensively. Additionally, changing the constraints with respect to computer time and property evaluations will likely lead to different results calling for further efforts towards better benchmark standardization and the introduction of various benchmark scenarios. We believe that explicitly considering the number of evaluations and computation times are still important to make the use of molecular design algorithms more accessible and relevant to the scientific community widespread application. Nevertheless, we think that useful insights are still to be garnered from our preliminary benchmark results. When looking at all the outcomes combined, one curious observation is that, unlike many previous benchmarking efforts, we do not find one algorithm to perform the best, or at least close the best, across all the tasks considered. This suggests that the newly introduced design objectives differ significantly in their structure requirements and algorithmic demands. Thus, we believe that different domains of chemistry and material sciences have different structure and design requirements that could potentially motivate the development of algorithms tailored to specific problem domains. At the very least, it suggests that significant further developments are necessary in the field in order to create a true champion algorithm for the design tasks considered.
long times detrimental to
Next, we will investigate the results of each of the molecular design sets in detail. The benchmark tasks inspired by the design of OPVs prove extremely difficult as most models fail improve upon the reference structures and, if they do, the improvements are very small. This suggests that the provided dataset already contains almost optimal solutions with respect to the fitness functions chosen. Nevertheless, the models should at least be able to reproduce the dataset results. We suspect that these difficulties of some of the models is likely related to uncertainties in the learned structural space [15T]] Notably, the property distributions depicted in Figure S8 show that the multiobjective function mainly dominated by the PCE values. We believe that increasing the relative importance of the SAscore in the objective function could improve this benchmark task even further. Nevertheless, our results suggest that the design of OPVs based on the Scharber model and an explicit account of synthetic accessibility is delicate and proves challenging for many common molecular design algorithms. Based on the judgement of an expert chemist, the best structures proposed by the molecular design algorithms are likely all stable and synthesizable showing the
$11

---

importance of accounting for that explicitly in molecular design benchmarks. It is interesting to see that the two best structures found in each of the objectives are quite similar which is likely because they do not stray far from the best compounds already contained in the reference dataset. The best structures for the donor design task seem to prefer having a 2.1.a-lhia moiety while the acceptor design task prefers extended annulated r-systems.
Looking at the objectives based on the design of OLEDs, they are all practically relevant as they are directly based on the design of TADF emitters. The molecular design algorithms showing most consistent performance across the three tasks are JANUS and GB-GA. Interestingly, SELFIES- -AAE also shows good performance over all three tasks and outperforms SMILES- VAE very significantly which could suggest difficulties of the SMILES-VAL to learn the underlying property distribution of the training dataset. Notably, none of the molecular design algorithms used outperform the best molecule in the reference data for the oscillator strength benchmark. Additionally, only three of the models, namely SELFIES VAE, JANUS, and GB-GA, deliver better properties for the combined objective than the reference structures possess. This demonstrates the difficulty of the corresponding benchmark tasks. Looking at the best performing molecules proposed by the models, the best structures are largely stable and synthesizable. However, there are a few features that are not ideal. In particular, the two molecules with the highest oscillator strength have an oxidized benzene ring that is likely reactive. Additionally, the molecule with the highest combined objective is present in its enol tautomer but it is likely that the corresponding ketone is preferred. This suggests that incorporating the SAscore a constraint helps increase the stability and synthesizability of the generated molecules but also suggests that additional filters are necessary for higher robustness in that regard. Moreover, some structural features required for favorable properties can be derived. In particular, for high oscillator strength values, extended x r-systems are beneficial which is well-known in the field of organic electronic materials [52 153].
Next, we need to inspect the outcomes of the protein ligand design benchmarks, which also have the lowest number of specific objectives, and is currently the only benchmark set without a multiobjective molecular design task. Most interestingly, while the general setup of all the three objectives is identical, we find the performance of the molccular design algorithms to be varying significantly depending the specific target. This suggests that designing ligands based molecular docking for the three proteins chosen differs significantly in terms of structural requirements. Most models succeed at proposing ligands with better docking scores than the best reference compounds for all three target proteins. In terms of the constraints, surprisingly, the VAES show the poorest performance which hints towards difficulties during training. Otherwise, JANUS also generally shows somewhat lower success rates which is likely caused by its limited features to bias structure generation by reference compounds and its tendency to explore more diverse structural regimes to potentially find new candidate solutions. This is relevant as the success rates are not only measured on the west-performing molecules but over all the structures proposed. In contrast, the LSTM-HC models, REINVENT and GB-GA consistently reach higher success rates demonstrating that they succeed at learning to propose desirable structures based on the reference dataset. Additionally, it likely also implies that exploration far outside of the distribution of the reference dataset is avoided. Thus, overall, REINVENT and GB-GA provide the most promising results in these protein ligand design benchmarks. When looking at the best-performing structures, it is notable that molecules resembling steroids seem to be good binders for at least two of the protein targets selected. Notably, there is one important aspect to be discussed that relates specifically to the dataset of reference compounds and the simulation workflow. While the DTP Open Compound Collection is particularly attractive to be used as training dataset because all these structures have been synthesized and tested experimentally for potential use as drug candidates, it suffers a major drawback in terms of data quality. Unfortunately, the majority of structures lack proper definition of stereochemistry despite often having multiple stereogenic centers. While our computational workflow is deterministic with respect to the particular stereoisomer produced when provided with SMILES without defined stereochemistry, it is still very likely that stereoisomers were produced and simulated that do not correspond to the ones of the actual molecules in the compound collection. Thus, when SMILES without well-defined stereochemistry are submitted to the property prediction workflow, the 3-dimensional structure generated in the process needs be inspected to derive the stereoisomer computed and guide subsequent experimental validation. For even more realistic drug design benchmarks tasks in the future, a reference dataset of structures that have both been experimentally synthesized before and
$12

---

have well-defined stereochemistry that reflects the experimental samples needs to be used.
When investigating the results of the chemical reaction substrate design benchmark tasks, we find that only the two GAs employed consistently outperform the best reference structures. This is particularly notable as both JANUS and GB-GA achieve large improvements over the best properties in the training set for all but the fourth benchmark objective showing that there is, in principle, significant room for further improvement. Additionally, in many cases, the proposed molecules perform significantly worse than the best- berf-penf reference compounds suggesting that many of the algorithms used seem to have difficulties to learn from the provided dataset. This could potentially point towards peculiarities of the provided training set as was the only one that was created in random way from the necessary core structure. It might be too small in size paired with a high diversity of structural moieties making it challenging to extract prrppere-property relationships. Notably, as expected because it goes against the Sil---van--PPoayy principle, the fourth benchmark objective proves to be extremely difficult with the best results providing only comparably small improvements over the reference set. The structures of the best proposed compounds for each of the design objectives provide some preliminary insights into the structural features effecting favorable properties (cf. Figure S4p. The sites donating hydrogen atoms seem to prefer being more substituted for both lowering the barriers and lowering the reaction energies. This effect can be understood by reduction of steric repulsion between these substituents and the bridge substituents when going from tetrahedral to trigonal geometry around the respective carbon atom. Particularly effective for lowering reaction energies seem moieties suitable for conjugation with the newly generated double bond in the product, owing to the stabilizing effect of extended conjugation. To break the correlation between reaction energies and activation energies, the introduction of two mesomerically electron- alectron-donating groups appear most suitable, potentially due to the opposing factors of sterics, lowering the reaction energy as outlined above, and electronics. Notably, deciphering substituent effects for these type of dyotropic reactions is known to be difficult [154), and would require more extensive and systematic analysis and corroboration which is outside the scope of this work.
The final benchmark of TARTARUS to be discussed is the model timing comparison. This is important as it can provide information about whether improved optimization performances originate from the utilization of more computational resources and longer training times. It should be noted that the absolute computation times reported are not useful as they strongly depend on the computer hardware available. Hence, we recommend them to be used relative rather than an absolute measure. Thus, when new models are being tested for timing, some of the molecular design algorithms provided here as reference should be run again with the same hardware as internal standards permitting relative timing comparison to all the other algorithms. Looking at the timing results, apart from the VAEs, we consider training to be still affordable for the deep generative models tested, even without the availability of GPUs. However, without GPUs, training the VAEs is computationally too expensive and detrimental to its widespread application, and even for the LSTM- HC models is comparably cumbersome. Additionally, while molecule sampling time, taking roughly on the order of a few minutes in most cases, is generally reasonable, we believe that there is considerable room for further improvement. This is particularly evident from the extremely good timing performance of REINVENT taking less than one minute for that step. Nevertheless, as it more important to propose several meaningful structures than to generate large number of structures fast, we believe that the sampling times of the all the algorithms tested are sufficient for practical purposes.
S5
Computational Details
S5.1
Molecular Design Algorithms
S5.1.1
VAE
We implemented variational autoencoder (VAE) architecture from the literature (31. In particular, the encoder used three ID convolutional layers with filter sizes of 9, 9, and 10 convolutionn kernels, followed by one fully connected layer that encodes a latent space of size 292. The decoder fed into three layers of gated recurrent unit (GRU) networks [107] with a hidden dimension of 501 each. Based on the respective training dataset to be used and the encoding used for representing
$13

---

the molecules (i.e., SMILES or SELFIES), we modified the alphabet size for converting the list of molecules into one-hot encodings. Additionally, the string length of the largest molecule from the dataset is used to pad sequences to the same length. Based on the dataset, these two parameters lead to minor modifications in the sizes of the first and last layers of the encoder and decoder, respectively. For each dataset, the models are trained using 24 CPUs and single GPU, with training times ranging between 2 and 10 hours. Generally, 80% of the reference molecules are used for training, while the remaining ones are used for testing. We manually optimize the learning rate, learning rate decay rate, number of epochs for training, and the batch size based on performance on the test set. After obtaining fully trained model, we perform structure optimization. The best known molecule for given task, i.e., the one with the most favorable score, is converted into a one-hot hot encoding and passed through the encoder to the latent space. Subsequently. Gaussian noise is added to the latent space vector to produce a population of vectors. These vectors are passed through the decoder to produce molecules that are evaluated via the respective objective function. The best molecule resulting from this population is then used to repeat this procedure.
S5.1.2 JANUS
Throughout the benchmarks, we used the default implementation of JANUS thttps://githii com/ com/aspuru-guziin -uter-uugu---gggaAAUUe--
We used version 1.0.3 of SELFIES, and the default valence constraints implemented in that version of SELFIES were employed. Additionally, no structure filters were employed in the genetic operators. For initiation, the file sample_start.. smi.es................mm..eee..e txt was populated with the 1,000 best molecules from the corresponding benchmark reference dataset. The generation_sian sintiiinneeeee parameters and the other parameters that influence the generations are modified based on the respective benchmark task. An artificial neural network (ANN) classifier is used for additional selection pressure for all the tasks (cf. JANUS+C in the original publication
S5.1.3 GB-GA
We employed of the implementation provided by Jensen 34 in the following repository: https: //github. thub. jttp//jjnnnnnuuppe CatggrrrppGGBBCCaa CB| For initiation, the file ZINC_first_In firrs_11000..t_ smi was pop- ulated with the 1,000 best molecules from the corresponding benchmark reference dataset for pre- conditioning to derive both mutation and crossover probabilities. The file scoring scoring_functions functions.....m py was modified each time to incorporate the corresponding fitness functions from different benchmark tasks.
S5.1.4 LSTM-HC
The long short-term memory hill climbing (LSTM-HC) approach uses pre-trained long short-term memory (LSTM) (105) recurrent neural network (RNN) to generate sequences of molecular strings 33]. and can be implemented with either SELFIES or SMILES. For each molecule, the string representation characters are transformed into one-hot encoded vectors, which are padded up to the length of the longest string in the respective reference dataset for each benchmark task. Similar to how it was implemented in the literature [33), the model consisted of 3 stacked LSTMs, each with 1024 hidden dimensions and dropout ratio of 0.2, followed by linear layer that outputs the logits of the next predicted character in the string sequence. That way, the output returns the probability for observing certain subsequent character, given the preceding characters in the sequence. The model was pre-trained with a batch size of 128, using the ADAM optimizer [1555 with a learning rate of 0.001. For optimization in each of the benchmark tasks, the top 2 known molecules with the best fitness values are used to create seeds for the LSTM structure generation. For SMILES, the strings are truncated randomly 25% to 75% of the size of the original sequence. For SELFIES, the strings are first randomized 5 times using reordered SMILES that represent the same molecule before they are being truncated. Subsequently, the truncated strings are passed into the model and the next characters are generated by sampling from the probability distribution that is generated by the LSTM output. Each iteration, the model generates 500 new molecules for evaluation, and the top 2 molecules are selected for the subsequent iteration. Additionally, the model is re- trained after each iteration on the set of new molecules for 10 epochs at a lower learning rate of 0.0001.
S14

---

S5.1.5 MoFlow
We employed the MoFlow framework developed by Zang et al. [32 for our molecular structure generation tasks. The implementation we used is publicly accessible https://githut. com/ We trained the model on 80% of the respective dataset, using batch size of 256 for maximum of 200 epochs. All model parameters were selected based on the example provided in the GitHub repository. Specifically, the B-Block of the model utilized 10 flows and had 512 hidden channels both layers. the A-Block, the model was configured with 38 flows and had 256 hidden graph neural network channels, along with hidden linear layers of sizes 512 and 64, Masking parameters were set with both a row size list and a row stride list of 1. A noise scale of 0.6 employed to introduce stochastic element into the behavior of the model. To optimize properties, we implemented hill-climbing-bart was randomly sampled, and the top 250 molecules were selected from total of 500 oracle calls for subsequent training. This cycle was repeated for up to 10 iterations, depending on the benchmarking task.
calvin-zcx/mof
algorithm. In this approach, the model
S5.1.6 REINVENT
In REINVENT 30), similar to LSTM-HC, an LSTM that is pretrained on a reference dataset is used prior for generating new SMILES sequences. Reinforcement learning (RL) is used to optimize for a particular fitness value. The LSTM model acts as an agent, where the addition of a character is an action, and the output probability distribution is the action probability, or the policy. The action probability is augmented by a scoring function S(m) E [-1, 1]. with larger scores corresponding to more desirable fitness values. SMILES that do not correspond valid molecules are assigned score of The algorithm updates the policy to increase the expected fitness while still keeping the learned conditional sequence probabilities produce valid SMILES that resemble the strucutural distribution of the reference dataset. The scoring function is a modified sigmoid function parameterized by the known fitness values from previous iterations. Let the set of known fitness values be The scoring function for molecule m is defined as
2
S(m) = + e-M(Sa) 1,
where a is the average of F, and b is the slope of the sigmoid. This sigmoid is defined as
1 In
2 cc+1 (021-1) -1)
b max(FF)
where e a threshold value set to 0.8. This means for the current fitness values in that the best known fitness will map to 0.8 in the scoring function. The smaller the threshold c, the stronger the reward for fitness values larger than the best one contained in F. Further details about the approach and model parameters are provided in the original paper [30].
85.2
Benchmarks
S5.2.1
Design of Organic Photovoltaics
Workflow To set up the property prediction workflow, we first implemented the Scharber model To do that, we constructed a simplified model to predict the short-cireuit current density (Jsc). The canonical way to do that is via definite integration of the spectrum of the light source over the energy range of the absorbed light of the simulated device [42]. To bypass having to perform the integration in the property simulation workflow and provide the Reference Air Mass 1.5 Spectra (AM1.5G) as part TARTARUS, we created regression model for the short-cireuit current density function of the band gap energy Eg. We used the AMI.5G spectrum that accounts for both direct and circumsolar irradiation and performed the integration via the trapezoidal rule as implemented in the numpy package. We found that simple Gaussian function of the form
41
Eg²
JsC A.e
(4)
S15

---

where A and B are fitting parameters, provides an excellent fit to the actual integration as demonstrated in the comparison between the original relationship derived from integration of the AM1.5G reference data (denoted "Original") and the fitted function (denoted "Fit") in Figure Importantly, the regression model assumes an external quantum efficiency of 0.65, which is typically used in the Scharber model [1] Typically, the band gap energy is approximated by the gap between the HOMO of the donor and the LUMO of the acceptor, which then, in conjunction with the regression model just described, allows to estimate the short-circuit current density that is needed to compute the power conversion efficiency (PCE) under the assumptions of the Scharber model [T
S6
Ongia
Figure S6: Short circuit current density as function of the band gap energy of the light-absorbing material in an organic solar cell. The original dependency results from the integration of the Reference Air Mass 1.5 Spectra (AM1.5G). The fit depicts the result of regression with a Gaussian function with two fitting parameters.
Besides the computation of the short-circuit current density, the procedure described in the original paper was followed [41]. The open circuit voltage (Voc) determined as energy difference between HOMO of the donor and the LUMO of the acceptor minus 0.3 eV of overpotential that is typically take place 2 Should the open circuit voltage be formally estimated to be negative then it is assumed to be O. In case the LUMO of the donor is less than 0.3 eV higher energy than the LUMO of the acceptor then the open eircuit voltage also assumed to be 0. The power conversion efficiency is estimated based on the equation
assumed be required for any charge separation
JSC P
PCE = 100% Voc FF.
(5)
where Vog is the open circuit voltage FF is the fill factor, Jsc is the short-circuit current density, and P... is the incident light intensity obtained by integration over the entire AM1.5G reference spectrum and therefore is constant. In the Scharber model, the fill factor is typically assumed to optimizable reach a value of 0.65 [156), which is thus inserted into this equation. For the first task of the organic photovoltaic design benchmarks, the goal is to design small organic donor molecule used with ww1----henll-66--buyyeee The LUMO energy level of PCBM is -4.3 eV and all light is assumed to be absorbed by the donor molecule. Notably, this is a common design objective that has been pursued the OPV literature [1 42 and was also implemented in the Harvard Clean Energy Project (CEP) [37, 40] For the second task, a small organic acceptor molecule to be designed that used with 7et.--cabbbzlleall-------- 70-di---thiieny-l-ee------------ 20,,,,0 10.30--ennothhadiazzoe (PCDTBT) as a donor. The HOMO energy level of PCDTBT is -5.5 eV and all light is assumed to be absorbed by the acceptor molecule. This benchmark task is based on a high-throughput virtual screening that was conducted in the literature to find new non non-fullerene acceptors (43)
acid methyl ester (PCBM) as acceptor.
The computational property prediction workflow accepts a SMILES string input and uses Open Babel [157] to generate initial guess Cartesian coordinates. Next, crest [188] is used for conformer search of the molecule via the iMTD-GC workflow at the GFN-FF level of theory [199] using the additional keywords -mquick, and -noreftopo. Subsequently, the lowest energy conformer resulting from this search used for geometry optimization at the GFN2-xTB level of theory (54). The properties of interest, i.e. HOMO energy, LUMO energy, HOMO-LUMO gap and molecular dipole moment, are obtained from the minimized structure the same level of theory.
S16

---

Importantly, to obtain better PCE estimates, we calibrated the HOMO and LUMO energies obtained from GFN2-xTB simulations against the corresponding density functional approximation (DFA) results taken from the Harvard Clean Energy Project Database (CEPDB) for the reference molecules of the CEPDB subset. In particular, the median HOMO and LUMO energies of all the DFA results contained in the CEPDB were used. Both the HOMO and the LUMO energies were calibrated against the respective medians the DFA HOMO and LUMO energies using a Theil-Sen estimator [61] as implemented in the python package eeiart- [162), due to its robustness with respect to outliers. The maximum number of subpopulations set to 250,000. The linear regression functions derived illustrated in Figure S7 and are provided in the following equations:
EioMMO EHOMO,GFN2- 0.8051 + 2.5377 eV
(6)
ÉLUMO, calsibrated M..UtOr,easbttteee
ELUMO,GFN2-2T2 0.8788 + 3.7913 eV
(7)
-2
JF
107
JAN (110-1
Figure S7: Two-dimensional histograms depicting the comparison of HOMO (left) and LUMO (right) energies obtained using DFA computations taken from the CEPDB against the corresponding energies at the GFN2-xTB level of theory obtained using the property simulation workflow developed for the design of organic photovoltaics benchmarks. The black lines were obtained via Theil-Sen estimators for linear regression.
As can be seen in Figure S7 the assumption of linear relationship between the energies at the GFN2-xTB level of theory and at the DFA level of theory is reasonable.
Dataset The CEPDB subset used for training all the generative models was obtained from the cwm/ Cww.///I/////a http/HIISS/neral-le ff waauuraalffiggge fingerprint, For all molecules in
following repository: this dataset, we simulated the corresponding properties of interest with the developed property simulation workflow (cf. Figure (8)).
S5.2.2 Design of Organic Emitters
Workflow The computational property prediction workflow accepts a SMILES string as input and uses both Open Babel [157] and RDKit [63 to generate initial guess Cartesian coordinates. Whichever guess coordinates have lower in energy the GFN-FF level of theory used as starting point for the subsequent conformer search of the molecule via the iMTD-GC workflow at the GFN-FF level of theory [159] using the additional keywords -mquick, ck, and -noreftopo. Subsequently, the lowest energy conformer resulting from this search is used for a geometry optimization at the GFNO- xTB level of theory (53 63 64.. Afterwards, the optimized geometry is used perform TD-DFT single point calculation at the B3LYP/6-31 IG level of theory (5) using the PySCF python package (version 1.7.6) 166) Density fitting utilized via the def2- aamversalllaie -jffii auxiliary basis set. For the excited state simulation, two excited states are generated for both the singlet and the triplet manifold. Notably, case the molecule contains iodine, the LANL2DZ basis set [167] is used for that atom. From the corresponding results, the energies of the lowest excited singlet and triplet states, and the oscillator strength for the transition between the ground state and the lowest excited singlet state are obtained. The singlet-triplet gap is derived energy difference between the first
S17

---

PCE1
SAS
PCE1 SSA
PCE1
PCE2
SAS
22 PCE2 SAS
PCE2
Figure S8: Two-dimensional histograms depicting the property distributions of the provided reference dataset obtained from the developed simulation workflow for the design of organic photovoltaics benchmark set. Top left: multiobjective target function PC PCEPCBM S Ascore against power conversion efficiency with PCBM acceptor PCEPCBM* top right: synthetic accessibility metric SAscore against power conversion efficiency with PCBM as acceptor PCEPCEM: bottom left: multiobjective target function PCE BT SAscore against power conversion efficiency with PCDTBT as donor PC PCEPcDI BT: bottom right: synthetic accessibility metric S'Accore against power conversion efficiency with PCDTBT as donor PCEPCDTBT.
10
Storet TrafLE
Si Siggks Irgie
Figure S9: Two-dimensional histograms depicting the property distributions of the provided reference dataset as obtained from the developed simulation workflow for the design of organic emitters benchmark set. Left: multiobjective target function +f12 AE(S)-Ti) - 1 against oscillator strength for the transition between S and So. fi2; middle: multiobjective target function +f12 -E(S-TT)-AAAE) (S 3.2 eV| against singlet-triplet gap, AE(S; - : oscillator strength for the transition between and S.. f12 against pogglet-tripler gap, AE(S TT).
3.2 eV|
S18

---

excited singlet and the first excited triplet state, Notably, only molecules with an SAscore below or equal 4.5 are subjected the property simulation workflow. All molecules above are assigned very unfavorable fitness.
originally comprising more than 970 million organic molecules, with over 400,000 structures possessing cycles and high degree of conjugation, that we herein refer as GDB-13sUB- These filters ensure that the molecules encompass cyclic T--sytems to increase representation of the relevant structural space for potential organic emitters. They were implemented via RDKit [163] and are summarized in Table S1O For the forbidden substructures, we utilized the following SMARTS patterns:
Dataset We developed comprehensive set of filters define subset of GDB-13 62
[C1, Br, I], ***c*, ***,
[c,c]N==, [0,, (0, "ta.:oo.....lla000oooa""""C"":""
*= *=[NH] *=N[[+;1R], [N,,,,,,,,,,s] (w...000..S.33311nn0000..... - [CH1] ++, *- [CH2] -*, *-[CH3]
Table S10: List of filters employed create the re-systems subset of GDB-13, GDB-13sUr, In each line, x denotes the value of the corresponding feature.
Number
Feature
Definition
Value
Charge
Charge of the molecule.
2
Radicals
Number of radical electrons.
3
Bridgeheud Atoms
Number bridgehead atoms.
4
Spiro Atoms
Number of spiro atoms.
0
5
Aromaticity Degree
Pereentage of aromatic non-hydrogen atoms.
0.5
6
Conjugation Degree
Percentage of conjugated bonds between non non-hydrogen atoms.
2: 0.7
7
Maximum Ring Size
Size largest ring.
8
Minimum Ring Size
Size of the smallest ring.
9
Suheructures
Presence of any forhidden substructures (see text).
FALSE
After filtering the GDB-13 dataset with these filters, and removing all molecules with an SAscore larger than 4.5, subjected the resulting molecules to the property simulation workflow described above to obtain the reference dataset for the design of organic emitters benchmarks. The corresponding property distributions are visualized in Figure
S5.2.3 Design Protein Ligands
Workflow To set up the computational property prediction workflow with docking calculations, we downloaded crystal structures for the three proteins of interest from the Protein Data Bank (PDB). The corresponding PDB codes are ISYH, 4LDE, and 6Y2F. The corresponding structures were all co-crystallized with ligand bound the protein. All these structures were processed using the protein preparation wizard [68] as implemented in the computer program Schrödinger (version 12.9). We created box around the bound ligand for each of the protein structures based re-docking accuracy. Subsequently, the native ligand was removed from the structure providing the starting setup for new ligands be bound in the property simulation workflow. This workflow accepts SMILES string as input, converts it into an initial 3D guess geometry using Open Babel [169] and exports it PDBQT file. Subsequently, the quality of the molecules checked by running obenergy. During reading of the SMILES, the corresponding structure is inspected for the presence of the following SMARTS patterns:
*1=**=*1, *1*=*=+1, *1**=*1, (F,C [Br]-C-C=OO,S,M,,
Cl, Br], [I], [S*x3] [SxX5]
[SæX6]
([F *=*=*=*, *= [NH]] [P,p] [F,, Cr], ss, CHC, C=C=C, [P.p] (=0)-++ =, C=C-N, NNN, "[*;R1]1~
000,
S19

---

109
Count
203
156
2000
253
25yh
20
Counts
117
060
4108 Scors
Figure S10: Two-dimensional histograms depicting the property distributions of the provided refer- ence dataset obtained from the developed simulation workflow for the design of protein ligands benchmark set. Top left: docking score to 4LDE A EALDE against docking score to 1SYH AESY H; top right: docking score to 6Y2F E against docking docking score to 6Y2F A 2 against docking score to 4LDE A FALLDE; bottom right: synthetic accessibility metric SAscore against quantitative estimate of wrug-likeness QED.
to ISYH A и: bottom left:
[#8] 1- [#6] 2 [#8] [#6] [#8] [#6] 12, N=C=0, C1CN1, [#6] (=#8]) [F, C1, I], [#6] (= [#8] - [#6] (- [[48]) [#6] (= [#8]) [#8], [(-[#6) -
These SMARTS patterns were developed from experience with the use of docking in conjunction with inverse molecular design algorithms. They correspond to reactive structural moieties that can lead to very favorable docking scores and should therefore be explicitly avoided 77) Additionally, it also tested whether the structure charged, possesses radical electrons, contains more than two bridgehead atoms, has least one ring with more than members and violates Lipinski's Rule of Five [128]. If at least one of these initial checks is true, the property simulation workflow is aborted and very bad fitness of 10,000 is returnd.. If all of these checks false, the workflow is continued. Subsequently. smina is used to add the ligand to the prepared protein structure of interest and perform docking simulation with an exhaustiveness setting of 10. The docked ligand structure with the lowest docking score is selected and the corresponding docking score returned by the simulation workflow. Following recent literature, 84 85] the quality of the generated docked structure was checked using obenergy and only values below 10,000 were accepted.
Dataset We obtained list of approximately 250,000 canonical SMILES from PubChem 170 that are part of the DTP Open Compound Collection 779, 80) Next, all the structures that do not pass any of the structural filters described in the previous section removed. Thus, obtained set of 152,296 molecules that are provided as reference dataset for the protein ligand design benchmark set and used it for training the generative models. For all molecules in this dataset, we simulated the corresponding properties of interest with the developed property simulation workflow (cf. Figure
S10
S5.2.4 Design of Chemical Reaction Substrates
Workflow The double hydrogen transfer in syn-sesquinorbornenes (cf. SigureSS11: was chosen as a suitable benchmark reaction due to its robust behavior in simulations. Additionaly, is well- described by the SEAM method as shown previously by Jensen [89]. For the benchmark, allow
$20

---

substitutions at the positions indicated by X or R in Figure S11B.
A
H
Figure S11: Intramolecular concerted double hydrogen transfer reaction syw-sesquinorbommenss (A) and positions in the syn-sesquinorbomene structure allowed to be substituted (B).
The general procedure of the reactivity workflow developed is outlined in the following steps:
1. Generation of initial conformation of the reactant.
(a) Try to generate coordinates with Open Babel [157].
(b) In case of failure, generate coordinates with RDKit [163
2. MMFF 171
optimization of the reactant.
3. crest [58] conformer search for the reactant via the iMTD-GC workflow.
4. Constrained MMFF optimization going from the reactant the product conformation consistent with reactant.
5. Full MMFF optimization the product.
6. crest conformer search for the product lowest conformation via the iMTD-GC workflow.
7. Crude reaction path interpolation with the geodesic interpolation method [72
8. SEAM 89 optimization of the TS.
9. Constrained crest conformer search of the TS via the iMTD-GC workflow.
10. SEAM re-optimization of the lowest energy TS conformer.
11. Calculation of reaction and activation energy from the SEAM model.
Prior to performing any property simulations, the molecule is checked for the presence of the syn-sesquinorbooreene core motif which is strictly required for the reaction of interest to occur:
[H] [cee] (*) [CQ;R2] ( [c00] [ce; (*) [C; (*)= [C;R1] (*) [C@;R2] (*) [[#;R2] [{:;R2]5) 5) [CEe]3 [+ 41] [Co] (+) 22] [c#0; R1] ([+]) [H]
Additionally, the input structure is inspected to ensure the absence of several reactive structural moieties described by the following SMARTS patterns:
[C-], [S-], [0-], [N-], '[++], [+-] [PH] [pн], [Nxx5],
*= [=,, [S8X3], [S*X4], [S8x5], [S8X6], [P.p],
[B,b,N,n,0,o,,,, [F, [I *=*=* ##* [0. [0..oS,,
",,,,..,,,s"",,,,,..........,, (00,o,8,,,N,;;IRR, *=N-(+;IR),
detailed description of some the steps in the property simulation workflow for the design of chemical reaction substrates is provided in the following bullet points:
1. SEAM optimization of the TS:
The SEAM model optimizes the transition state as the crossing point between reactant and product force fields. In this work, we used the GFN-FF force field [559, and shifted the product foree field reproduce the energy difference between the reactants and products at the GFN2-XTB level of theory [54). For the product force field shift, used the lowest
S21

---

12
-20
Figure S12: Illustration of the outlier detection cutoff for the elliptic envelope used for the property simulation workflow of the design of chemical reaction substrates benchmark set. Data points appearing outside the black line are regarded as outliers and assigned a very low fitness.
energy conformers of the reactant and product, respectively. The optimization was carried out on the upper SEAM surface by introducing a small coupling of 0.001 eV between the states. The guess structure for the TS was taken from geodesic interpolation between reactant and product geometries [72]. The TS itself was optimized with the geome TRCC package [73] via an interface to pyscf [66). The subsequent crest transition state conformer search used constrained C-H bond lengths with force constant of 1.0 a.u. in strength. All calculations were carried out using the POLANY] package, which will be reported in detail elsewhere. The activation energy was taken the difference of the GFN-FF energy of the optimized TS structure and the energy of the lowest energy conformer of the reactant. The reaction energy was taken the energy difference between the lowest energy conformers of the product and reactant, including the energy shift of the product force field. These energies are electronic energies from the force field. While free energies would be preferable, it would increase the complexity of the workflow, and significant error cancellation expected for this intramolecular reaction.
2. crest conformer searches:
We opted to perform the crest conformer searches with the keywords -mquick, gff2// Sfn2///fuff. ggfff, and -noreftopo. Although it would be more consistent use -gfnff, we found that it produced larger number of outliers in preliminary test calculations.
3. Outlier detection:
A Although the workflow is very robust, it does produce some outliers due to failures in the conformational search or the TS optimization. the worst cases, these outliers could display extreme properties such as artificially low or high reaction activation energies, and, hence, sometimes emerge as the best performering structures. One solution is increasing the exhaustiveness of the conformational sampling, but this also leads to significantly longer computational run-tima ain-time.. Instead, we opted for an outlier detection model based the computed reaction and activation energies. According to the Goll-Evans-Polanyy relationship [148), the reaction energy is expected to correlate with the activation energy. Therefore, outlier model based on the 21) scatter plot of these two properties can be from sccikit-learn set to 0.00035 by visual inspection to
used to identify potential outliers. We used the [17). The contamination hyperparameter allow for considerable deviation from the correlation observed in the reference dataset, while simultaneously capturing clear outliers. This outlier detection scheme illustrated in Figure ST2) Accordingly, we incorporated this outlier detection into the property simulation workflow and assigned an extremely low fitness of 10 to outlier molecules.
Dataset We constructed a dataset of 60,850 molecules by performing multiple mutations of various SELFIES representations of the parent syw-sesquinorborrenne. Specifically, the starting structure was used to produce 20 randomly reordered SMILES strings. The resulting SMILES were converted to SELFIES and 20 random mutations were conducted via CTONED-SELFIRE [12]. All mutants were subsequently converted to SMILES and checked for the presence of the core motif. Mutations that lead to structures without the core motif were discarded. Additionally, molecules that did not pass the filters were also removed. This procedure was repeated on all unique and valid mutated structures until the final number of molecules was reached. Notably, we used the default mutation settings of STONED-SELFIES 2T2). When running the dataset through the property simulation workflow, only
S22

---

35 molecules either did not run through the property simulation workflow properly or were outliers, corresponding to a success rate of 99.94%. The results of running the reference structures through the property simulation workflow with respect to both property distributions and computational run time are illustrated in Figure S13
2502
10003
Reaction Energy
Frequericy
3
130
Time (snc)
Activation Energy
Figure S13: Simulation results for the reference structures in the developed simulation workflow for the design of chemical reaction substrates benchmark set. Left: Histogram of the simulation workflow run times; right: two-dimensional histograms depicting the property distributions of the provided reference dataset: reaction energy AF, against activation energy AE
S6 Model Timing Comparison
All deep generative models (SMILES- VAE, SELFIES- VAE, MoFlow, REINVENT, SMILES-LSTM- HC and SELFIES-LSTM-HCC were trained using the first 80% of the molecules from the protein ligand design tasks. The remaining 20% of the dataset was used to assess convergence and perform manual hyperparameter optimization. The corresponding figure in the main text provides four timing benchmark metrics: (1) Training Time with GPU, (2) Sample time for 10,000 molecules. (3) Epoch time without GPU and (4) Epoch time with GPU. For each metric, 5 independent experiments were conducted obtain both average and standard deviations of the respective values. For the first metric. all the generative models were trained using 24 CPUs (AMD Rome 7532 2.40 GHz 256M cache L3) and single GPU (Tesla A100). For REINVENT, MoFlow, and the LSTM-HC models, a batch size of 512 was used for training. In contrast, the VAE models performed better in terms of validity and reconstruction with smaller batch size of 8. For sampling times of 10,000 molecules, we used fully trained deep learning models obtained after the training time measurements and determined the time required for generating 10,000 molecules in the presence of 24 CPUs (AMD Rome 7532 2.40 GHz 256M cache L3) and single GPU (Tesla A 110)). The batch size was 512 during sampling for all deep generative models. For GAs, the molecule sampling time was measured only in the presence of 24 CPUs as both JANUS and GB-GA lack GPU support. Notably, in the models that use SELFIES as a molecular representation, the time nceded for translating SELFIES to SMILES was included the final reported time. For JANUS and GB-GA, this metric obtained by increasing the population size to 10,000 and by using single generation runs and dummy fitness function that returns a value of for any molecule. For the single epoch times, we report the time required for training models in the presence of either 24 CPUs (AMD Rome 7532 @ 2.40 GHz 256M cache L3) or 24 CPUs (AMD Rome 7532 @ 2.40 GHz 256M cache L3) with GPU (Tesla A A100). For the deep generative models, a batch size of 512 was used. Notably, for both GB-GA and JANUS, the CPU epoch time corresponds to the total time the models need to precondition the genetic operators based on dataset of reference structures. Additionally, for both these models, preconditioning was performed by providing 1 1,000 random SMILES from the reference dataset of the protein ligand design task to the molecular design algorithms. Finally, using the CPU epoch times and the number of epochs needed for training when also using a GPU, we estimated the total CPU training times for all the deep generative models by multiplying these two numbers. Notably, we decided not to estimate the corresponding estimated errors via error propagation as these values are merely rough estimates.
$23


Reviews:
[{"review_id": "9GwL8ieXEK", "reviewer": "Reviewer_uBzk", "summary": "", "questions": "", "limitations": "I don’t see much effort to rationalize what could have been the origin of performance differences between different generation tasks. While these models (VAE) are inherently complex, their relative strengths and weakness should be commented to get a better feel of the trends reported through different tables in the manuscript.", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "In a time when paradox of choice is gripping every domain, what kind of deep generative model may become suitable for a given task is an important problem. Authors have competently executed their to address this important problem and came to the conclusion that the degree of sophistication of the deep generative model (or the manner in which chemical featurization/encoding is carried out) does not necessarily guarantee stronger performance in inverse design problems.", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": ""}, {"review_id": "MlCHCAMrfu", "reviewer": "Reviewer_zwfp", "summary": "", "questions": "", "limitations": "The authors didn’t adequately discuss the negative social impact of their work. For example, the molecular generative models might be used to generate chemical weapons and hazardous materials.", "rating": 7, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "- The paper develops a simulation-based workflow to evaluate generated molecules, employing openbabel, Crest, XTB, etc. These evaluation workflows are much more robust than widely used LogP, synthesizability scores, etc. \n- The benchmark tasks cover a broad range of scientific domains, including organic photovoltaic design, protein ligand design, chemical reaction subtract design. \n- It also extensively benchmarked 7 different molecular generative model algorithm\n", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": ""}, {"review_id": "Dai0Hoxc8d", "reviewer": "Reviewer_yzPs", "summary": "", "questions": "", "limitations": "Yes.", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "The proposed evaluation methods based on physical simulation are all sound. And the evaluation is aligned of the real-world application.  \nThe analysis of some selected molecule design algorithms reflects that there is still a lot of root for developing more useful methods.  \nThe codebase of this benchmark is well-developed, which provides a good standard to follow up for researchers in this community.", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": ""}, {"review_id": "rl6wBXPXQl", "reviewer": "Reviewer_8Sdd", "summary": "", "questions": "", "limitations": "Yes.", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "(1) The presented TARTARUS provides four interesting and useful datasets for developing inverse molecule design methods. The tasks of these datasets are tightly related to realistic and practical problems in real-world applications. These datasets can be broadly useful for researchers in the field of inverse molecule design.  \n(2) The writing of this paper is clear and well-organized. The details of the creating process of all datasets are clearly presented.", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": ""}]

Author rebuttal:

Meta-review: ## Summary

The authors curate training sets and develop simulation-based evaluation suites for several interesting generative molecular design problems. They then evaluate a set of generative models. 

## Strengths

- [quality] The codebase of this benchmark is well-developed, which provides a good standard to follow up for researchers in this community.
- [significance] The analysis of selected molecule design algorithms reflects that there is still a lot of root for developing more useful methods.
- [quality] The proposed evaluation methods based on physical simulations are all sound. And the evaluation is aligned of the real-world application.
- [significance] The benchmark tasks cover a broad range of important scientific domains, including organic photovoltaic design, protein ligand design, chemical reaction substrate design.
- [clarity] The paper is very easy to read. 

## Weaknesses

- [quality] Not comparing to experimental values weakens the conclusions
Decision: Accept (Poster)