# Repurposing quantum chemical descriptor datasets for on-the-fly generation of informative reaction representations: application to hydrogen atom transfer reactions

Javier E Alfonso-Ramos, Rebecca M Neeser, Thijs Stuyver

► To cite this version:

Javier E Alfonso-Ramos, Rebecca M Neeser, Thijs Stuyver. Repurposing quantum chemical descriptor datasets for on-the-fly generation of informative reaction representations: application to hydrogen atom transfer reactions. Digital Discovery, 2024, 3 (5), pp.919-931. 10.1039/D4DD00043A. hal-04951875

HAL Id: hal-04951875

https://hal.science/hal-04951875v1

Submitted on 17 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Cite this: DOI: 00.0000/xxxxxxxxxx

# Repurposing Quantum Chemical Descriptor Datasets for on-the-Fly Generation of Informative Reaction Representations: Application to Hydrogen Atom Transfer Reactions

Javier E. Alfonso-Ramos,ᵃ Rebecca M. Neeser,ᵇ,ᶜ and Thijs Stuyver*ᵃ

Received Date
Accepted Date

DOI: 00.0000/xxxxxxxxxx

In this work, we explore how existing datasets of quantum chemical properties can be repurposed to build data-efficient downstream machine learning models, with a particular focus on predicting the activation energy of hydrogen atom transfer (HAT) reactions. Starting from a valence bond (VB) analysis of a generic HAT process, a set of informative descriptors is identified. Next, a surrogate neural network model is constructed to predict an informative representation, based on the identified VB descriptors, with the help of a publicly available dataset of (pre-computed) quantum chemical properties of organic radicals. We demonstrate that coupling the resulting on-the-fly informative representation to a secondary machine-learning model for activation energy prediction outperforms various predictive model architectures starting from conventional machine-learning inputs by a wide margin, at no additional computational cost. As a bonus, by basing their final predictions on physically meaningful descriptors, our models become inherently interpretable. Finally, because of the extreme data efficiency of our descriptor-augmented models, we are able to fine-tune and apply them to small datasets across various reaction conditions, settings and application domains, ranging from regular (liquid phase) synthesis, over metabolism and drug design, to atmospheric chemistry.

## 1 Introduction

In recent years, the interest in predictive (machine learning) models for chemical reactivity has soared. For tasks for which hundreds of thousands to millions of data points are available, e.g., forward and retrosynthesis prediction, accurate and robust machine learning (ML) models have been designed. In the absence of data scarcity, neural network architectures that learn from simple structural descriptions of molecules can be adapted to these respective tasks with great success. For example, both graph neural networks (GNNs) operating on 2D molecular graphs and Transformer models operating on Simplified Molecular Input Line

Entry Specification (SMILES) ¹ strings can be trained to predict the main product of chemical reactions with roughly 90% accuracy on the popular USPTO benchmark dataset.²⁻⁴

For more specialized reactivity tasks (the prediction of activation energies/reaction rates, yields, enantiomeric excess etc.), however, data tend to be scarce, and often number in the hundreds rather than in the hundreds of thousands. A straightforward application of the same type of data-hungry model architectures is generally not feasible in these cases, as training these models from scratch means that they must learn how to properly interpret a molecular structure. One potential approach to improve the data efficiency of machine learning models is introducing domain knowledge about the underlying physics/quantum chemistry into the model. Models starting from physically meaningful input representations, i.e., representations based on quantum mechanical (QM) reactivity descriptors, have been demonstrated to reach reasonable accuracy for datasets consisting of only a couple hundred or a few thousand data points,⁵⁻¹⁴ and they have been reported to be significantly more generalizable, i.e., they are significantly more accurate in out-of-sample predictions, than conventional, structure/graph-based, analogs.⁷,⁸,¹⁵

Unfortunately, QM descriptors tend to be computationally ex-

<footnote>ᵃ Ecole Nationale Supérieure de Chimie de Paris, Université PSL, CNRS, Institute of Chemistry for Life and Health Sciences, 75 005 Paris, France; E-mail: thijs.stuyver@chimieparistech.psl.eu</footnote>

<footnote>ᵇ Massachusetts Institute of Technology, Department of Chemical Engineering, 02139 Cambridge (MA), United States.</footnote>

<footnote>ᶜ ETH Zürich, Institute of Pharmaceutical Sciences, 8093 Zürich, Switzerland.</footnote>

<footnote>† Electronic Supplementary Information (ESI) available: [details of any supplementary information available should be included here]. See DOI: 00.0000/0000000.</footnote>

<footnote>‡ Additional footnotes to the title and authors can be included e.g. 'Present address:' or 'These authors contributed equally to this work' as above using the symbols: ‡, §, and ¶. Please place the appropriate symbol next to the author's name and include a \footnotetext entry in the correct place in the list.</footnote>
pensive – since they typically involve an elaborate workflow con-
sisting of conformer generation, geometry optimization and fi-
nally a single-point density functional theory (DFT) calculation
for each unique molecule – limiting their applicability in broad
virtual screening campaigns. A promising strategy to circum-
vent this limitation is to use separate machine learning mod-
els to predict these QM reactivity descriptors on-the-fly, start-
ing from a simple 2D (graph) or geometry-based representation.
Some prominent examples of such descriptor prediction models
are Chemprop, DeepMoleNet, SchNet and AIMNet<sup>16–19</sup>. In 2020,
Guan et al. proposed to directly integrate such a descriptor pre-
diction model into a reactivity prediction model as a surrogate.<sup>20</sup>
In their approach, a directed message-passing neural network (D-
MPNN) model<sup>16</sup> was trained to predict various atom- and bond-
level QM descriptors from input molecular structures. These QM
descriptors were subsequently used to augment the hidden repre-
sentations of a Weisfeiler-Lehman (WL) based GNN<sup>21</sup> to predict
regiochemical preferences for (aromatic) substitution reactions.
For training sets as small as 200 data points, the resulting "QM-
augmented" models reached accuracies of 85% and more on a
random data split, outperforming regular GNN analogs by more
than 10%. Furthermore, the developed models were extremely
fast; individual predictions required 100 ms or less, implying
that they could be readily employed in high-throughput screen-
ing campaigns. In a follow-up study, Stuyver et al. adapted this
model to predict activation energies for competing (gas-phase)
E2 and S<sub>N</sub>2 reactions and observed a remarkable improvement
in accuracy for their GNN model upon QM-augmentation as well,
particularly in the data-limited regime.<sup>22</sup>

Even when a surrogate model strategy is viable, initial training data for the model still needs to be generated, and in some situations, this may actually be more expensive than the explicit calculation of QM descriptors at inference time. For example, in their proof of principle study, Guan et al.²⁰ generated a dataset of selected QM properties for 130k organic molecules, extracted from the ChEMBL²³ and Pistachio²⁴ databases.

A potential strategy to circumvent the expensive QM descrip-
tors calculations altogether would be to use pre-existing datasets.
In recent years, various general-purpose public datasets of com-
puted QM properties of molecules have been released.²⁵⁻³¹ While
these datasets have found dual use as benchmarking targets,³²⁻³⁴
as well as training data for generative models,³⁵⁻³⁷ their repur-
posing for (indirectly related) downstream predictive tasks, e.g.,
the prediction of chemical reactivity, has remained largely unex-
plored up to this point.

Here, we will demonstrate that such a data repurposing strat-
egy can be readily applied to efficiently learn activation ener-
gies for hydrogen atom transfer (HAT) reactions. Starting from a
generic description of HAT reactivity within the Valence Bond per-
spective, we determine a set of qualitative characteristics that can
be expected to capture the main trends for this type of reactivity.
Based on this analysis, we identify the recently released dataset
of QM properties for 40,000 closed-shell organic molecules and
200,000 radical analogs by Paton and co-workers, BDE-db,²⁹ as
a useful source of informative descriptors and build a surrogate
model able to infer these descriptors on the fly across a broad

chemical reaction space based on a simple SMILES string input.²⁹
The QM representation emerging from this surrogate model is
subsequently demonstrated to outperform traditional structure-
based representations for various machine learning model archi-
tectures on a challenging, in-house generated dataset of diverse
hydrogen abstraction reactions by a large margin. Most impor-
tantly, models based on the developed surrogate QM represen-
tation are extremely data efficient, approaching convergence for
training dataset sizes of only a couple hundred points.

Finally, we demonstrate the practical usefulness of the devel-
oped representation by considering various real-world applica-
tions. HAT reactivity plays a central role in various sub-fields
of chemistry. For example, HAT reactions constitute the starting
point of many oxidation pathways in polymer, atmospheric and
combustion chemistry. ³⁸⁻⁴³ Additionally, they are key transforma-
tions giving rise to the breakdown of active pharmaceutical ingre-
dients, so their kinetic and thermodynamic properties facilitate
the prediction of the stability of these compounds, ⁴⁴,⁴⁵ as well
as the formation of potentially hazardous metabolites. ⁴⁶ Finally,
HAT reactions also play a central role in various common synthe-
sis strategies, e.g. in several flavors of C-H functionalization and
photocatalysis, and a lot of effort has been devoted to gaining in-
sights into the factors that affect their rates and yields. ⁴⁷⁻⁵⁰ As a
result, several small datasets have been produced for this reaction
class in recent years, each with distinct reaction scopes and con-
ditions, depending on the specific application for which they were
designed. ³⁸,⁴⁶,⁵⁰⁻⁵⁷ Because of their heterogeneity, these datasets
cannot be combined easily. As will be demonstrated below, with
our extremely data-efficient model architectures, we are able to
meaningfully learn from these small datasets, and generate ac-
curate predictions across the chemical spaces associated with the
respective application domains, at both an unprecedented speed
and an extremely low computational cost.

In essence, the proposed strategy rests on making the over-
all machine learning model architecture modular, with a core
informative representation prediction part, followed by a fairly
linear and extremely data-efficient downstream reactivity model
part, that can be easily retrained on/finetuned to the specific
application or dataset at hand. Overall, we believe that the
presented work provides a blueprint for how – starting from
physics/chemistry-inspired mechanistic insights – useful (exist-
ing) datasets of QM properties can be identified and leveraged for
downstream predictive tasks. Furthermore, it underscores the im-
portance of community efforts to generate and curate QM prop-
erty datasets – and to include as many informative descriptors as
possible in them. At the end of this contribution, we briefly reflect
on the transferability of the approach to other reaction types.

## 2 The Valence Bond reactivity model as a source of inspiration for descriptor selection

As indicated in the introduction, to gain qualitative insights into
the main factors driving HAT reactivity, we took inspiration from
qualitative Valence Bond (VB) theory.⁵⁸,⁵⁹ In a VB reactivity anal-
ysis, reactivity diagrams, reflecting the shape of the potential en-
ergy surface (PES) associated with a given chemical system along
a chosen reaction coordinate, are constructed. Usually, one starts by examining the evolution of the so-called “diabatic curves”, which represent the wave function associated with a limited number of VB structures; i.e., they correspond to a localized VB configuration (e.g., a specific Lewis structure of the reactants or products). These individual diabatic curves interact along the reaction coordinate, and collectively they give rise to the “adiabatic” curve, i.e., the full ground-state PES associated with the chemical system under consideration. The transition state then corresponds to the top of the hill of the adiabatic state, which usually lies in the same region as the crossing point between the reactant and product diabatic states. An in-depth discussion of the theory behind the construction VB reactivity diagrams is provided in Ref.⁵⁹. In Figure 1, a simplified sketch of such a diagram is presented for the model reaction between H₂ and H·.

<figure><img src="image_2.png" /><figcaption>Fig. 1 An example VB reactivity diagram for a model HAT reaction between H₂ and H·. Along the reaction coordinate, the reactant diabatic state (blue) rises in energy (since the hydrogen atoms contributing electrons to the formal covalent bond on the left become increasingly separated as the reaction proceeds), the product state (red) drops in energy, since the product covalent bond becomes increasingly stabilized as the associated hydrogens approach. The adiabatic state (black) is the full ground- (or adiabatic) state resulting from the mixing of the localized diabatic states.</figcaption></figure>

From a generic Valence Bond diagram, the following universally valid – yet approximate – expression for the barrier height can be derived:⁵⁹

$$ \Delta E = f * (G_r + G_p) + \frac{1}{2} \Delta E_{rp} + B \quad (1) $$

where $G_r$ is the promotion energy on the reactant side, i.e., the energy difference between the reactant and product diabatic states in the reactant geometry, $G_p$ is the corresponding promotion energy on the product side, $f$ depends on the curvature of the diabatic states and is a fraction, i.e., a number between 0 and 1, and $\Delta E_{rp}$ corresponds to the thermodynamic driving force, i.e., the energy difference between reactants and products. These first

two terms in the expression collectively indicate the height of the crossing point between the reactant and product diabatic states. The third and final term, $B$, is a correction term corresponding to the interaction between the diabatic states in the transition state geometry, i.e., it corresponds to the energy lowering from the crossing point of the reactant and product diabatic states and the top of the hill on the adiabatic state.

Eq. 1 is an essentially linear – and inherently interpretable – one, and has been applied successfully to various areas of chemical reaction space, leading to robust and fairly accurate barrier estimates.⁶⁰⁻⁶⁷ Unfortunately, obtaining reliable values for each of the featured parameters is ambiguous and computationally expensive for some reaction types. For HAT reactions, however, most of the terms in the expression can straightforwardly be connected to unambiguous QM descriptors, which can be readily computed with high accuracy, hence our decision to focus on this type of reactivity in this proof of concept study.

More specifically, in HAT reactions, the promotion energies, $G_r$ and $G_p$, are proportional to the vertical/frozen bond dissociation energy (BDE) – or, more accurately, the bond dissociation free energy (BDFE) – of the covalent bond being broken in the reactant and product geometries respectively.⁵⁹ The thermodynamic driving force, $\Delta E_{rp}$, in its turn corresponds to the (relaxed) BD(F)E by definition. $f$ and $B$ have been demonstrated to be fairly constant within narrow reaction families, e.g., when considering a carbon-based radical abstracting H atoms from a C-sp³ bond.⁵⁹ Since our ambition here is to generalize across the full space of (neutral) C, H, N, O containing organic radical reactions, this assumption can be expected not to hold perfectly. It has been demonstrated however that both of these parameters are connected to a variety of descriptors that can in fact be quantified straightforwardly.

First and foremost, $f$ depends on the electrostatic interactions along the *abstraction site - hydrogen - radical site* point-charge array in the transition state (TS) geometry, which can be probed indirectly by considering atomic charges (in the reactant and product geometries).⁶⁸ Additionally, $f$ also depends on the sterics around the radical and abstraction sites – steric bulk around those sites can be expected to strain the transition state geometry and hence should result in a higher crossing-point between the reactant and product diabatic curves. Commonly, buried volume ($V_{bur}$) is used as a descriptor to quantify steric bulk around reactive sites.⁶⁹ Finally, $f$ has also been suggested to depend on the extent of delocalization in reactant and product, which can be probed through consideration of the spin density ($\rho$).⁵⁹ However, an alternative viewpoint is that losses and/or gains in delocalization energy may already be fully reflected in the thermodynamic driving force term ($\Delta E_{rp}$), so that the impact of delocalization on $f$ may actually be negligible.⁶⁴ As indicated above, $B$ depends on the interaction between the various diabatic states, and has been connected primarily to the (local) ionization potentials and electron affinities – and/or electron negativities – at the radical/abstraction sites.⁵⁹

The discussion in the previous paragraphs provides us with an ideal QM descriptor representation that should render the mapping function from reaction to HAT barriers as linear as possible. As already indicated in the introduction, computing this complete
input representation would be complex and expensive – many dif-
ferent DFT calculations would be needed (for both charged and
uncharged species) for every data point. At the same time, the lit-
erature suggests that even a subset of the descriptors mentioned
above could be sufficient to build an informative representation:
for small subsets of the chemical space we are interested in, vari-
ous groups have previously constructed effective multivariate lin-
ear models, as well as some more advanced architectures (ran-
dom forests, neural networks, etc.) with a subset of the descrip-
tors listed above. 46,51,52,70

As already indicated above, encouraged by these findings, we
selected the BDE-db dataset as a suitable source of descriptors.
This dataset contains (relaxed) BDFEs – at room temperature
in the gasphase – and BDEs, (Mulliken) atomic charges, and
(Mulliken) spin densities for over 200,000 organic radical com-
pounds.²⁹ Additionally, since accurate geometries are included
in this dataset, we were able to expand this descriptor list with
buried volume values at a marginal computational cost. Finally,
we also used these geometries to rapidly compute vertical, i.e.,
frozen or unrelaxed, BDE values at extended tight-binding (xTB)
level of theory (these final low-level descriptors were found to
be remarkably predictive – even at this low theory level – for ac-
tivation energies of HAT reactions between P450 enzymes and
organic compounds in a recent study by Gingrich et al.⁴⁶ – see
the Methodology section for more details).

With this dataset-derived representation selected, we set out
to train a surrogate descriptor model, in combination with a top-
level reactivity model, to probe HAT reactivity on the fly.

3 Methodology

**3.1 Extracting and curating the descriptor data from the**
**BDE-db dataset**

As indicated in the introduction, the training data for the surro-
gate neural network (vide infra) was extracted from a dataset re-
cently released by Paton and co-workers, consisting of QM proper-
ties for radicals and associated closed shell molecules of less than
ten heavy atoms, without formal charges on any of the atoms, and
containing exclusively C, H, N, and O atoms. We complemented
the registered descriptors, i.e., BDFEs, spin densities, and par-
tial charges, with Vbur values around the radical site and frozen
BDE values computed at GFN2-xTB level of theory. Vbur values
were computed with the help of the Morfeus package<sup>71</sup> with a
radius of 3.5 Å.<sup>72</sup> Frozen BDEs were computed by re-optimizing
the extracted geometries at GFN2-xTB level of theory and then
performing single-point calculations on the bond fragments, i.e.,
the organic radical and a hydrogen radical respectively. In to-
tal, the complete set of descriptors could be obtained for 186,240
dissociation reactions in this manner.

The entire dataset in its curated form can be accessed and
downloaded from https://figshare.com/articles/dataset/
dataset_surrogate_model_cs/24754341. It should be noted
that during the preparation of this manuscript, an expanded ver-
sion of the dataset, covering the elements S, Cl, F, P, Br and I as
well was released.⁷³ As such, the approach described here could
be expanded straightforwardly to a much broader scope of HAT

reactions in the near future.

3.2 Generation of the HAT dataset

To enable comparisons between various types of machine learn-
ing model architectures, a dataset of HAT reaction profiles was
generated by combining multiple bond dissociation reactions in-
cluded in the BDE-db dataset.

A fully automated reaction profile computation workflow based
on autodE⁷⁴ and Gaussian16⁷⁵ was set up to generate this
dataset, starting from an input list of reaction SMILES. This work-
flow is analogous to the workflow recently developed for cycload-
dition reactions by Stuyver et al.;⁷⁶ an in-depth discussion can be
found in Section S1 of the Supporting Information.

It should be noted that around 40 billion HAT reaction SMILES
can in principle be constructed by making combinations between
bond dissociation reactions in the BDE-db dataset. Most of
these reactions are very similar, involving HAT among two car-
bon atoms. We aimed to select a tiny sample (2,000 reactions)
from this huge search space, covering maximal structural diver-
sity, to compute explicitly. To this end, we iteratively sampled
enumerated reactions in such a way that the differential reaction
fingerprint (radius = 3; nbits = 2048)⁷⁷ distances with all previ-
ously sampled reactions exceeded 0.85 (see Section S1 for more
details).⁷⁸

As the level of theory for the DFT calculations, we selected
M06-2X/def2-SVP⁷⁹,⁸⁰ for optimizations and frequency calcula-
tions, i.e., thermal corrections. Single-point energy refinements
were performed at M06-2X/def2-TZVP⁷⁹,⁸¹ level of theory. The
choice for the functional and basis set was inspired by the bench-
marking results from Paton et al.⁸²

One technical aspect where we decided to deviate in our work-
flow from the previous work on cycloaddition reactions is the con-
former selection.⁷⁶ Previously, the default number of conform-
ers in autodE (300) would be generated for every species along
the profile, a root mean squared deviation (RMSD) threshold of
0.3 Å would then be applied to remove very similar conformer
geometries, after which the lowest energy conformer, computed
at GFN2-xTB level of theory, would be selected. Remarkably, for
HAT reactions, this methodology leads to poor reproducibility of
the activation energies (mean absolute errors between successive
runs on 30 randomly sampled reactions amounting to more than
2.1 kcal/mol for the activation energy, and ~ 0.5 kcal/mol for the
reaction energy; cf. Section S1). After some benchmarking, we
adapted these conformer selection settings, selecting 1000 con-
formers for every species, using an RMSD cut-off of 0.1 Å, and
ranking conformers based on single-point M06-2X/def2-SVP⁷⁹,⁸⁰
energies. In this manner, two successive runs of the workflow on
the same set of reaction SMILES result in an acceptable mean ab-
solute error (MAE) of ~1 kcal/mol for the activation energies and
~0.6 kcal/mol for the reaction energies.

Tunneling corrections were computed based on the semi-
classical approximation.⁸³ In this approximation, the potential
energy barrier is approximated by an Eckart potential, for which
the tunneling probability is computed analytically and the tun-
neling transmission coefficient is integrated numerically. More
details can be found in Section S1 and reference$^{83}$.

In total, 1511 reaction profiles were successfully computed in
this manner (75.6% success rate).

## 3.3 Surrogate neural network model for QM descriptor prediction

A graph neural network (GNN), consisting of a directed message-
passing neural network (D-MPNN) encoder and a feed-forward
neural network (FFNN), was selected for the fast prediction of the
QM descriptors. ¹⁶,²⁰ Inference for an individual molecule starts
by parsing its SMILES string, after which a graph-based input
representation is constructed, with atoms as nodes and bonds as
edges. The D-MPNN module subsequently encodes this molecu-
lar graph, resulting in learned atom and bond features. Finally,
the learned representation of the molecule is used as input to
the multitask FFNN, which yields the final predictions (cf. Figure
2). A more detailed description of the model architecture can be
found in Section S2 of the SI.

The surrogate model is trained to predict two atom-level de-
scriptors (atomic Mulliken charges and spin densities), as well
as three molecule-level descriptors (relaxed BDFEs, frozen BDEs,
and V<sub>bur</sub> of the radical sites). Based on a hyperparameter grid
search, the depth was set to 5 for the D-MPNN and 4 for the
FFNN. For the hidden size, 1200 neurons were selected for the D-
MPNN and 900 for the FFNN (cf. Section S2.1 in the Supporting
Information).

Once a reasonable set of hyperparameters was selected, the
model was trained on a new data split, where all the radical
species and molecules present in the constructed HAT reactiv-
ity dataset (vide supra) constitute the test set, and the remaining
points are split between training (80%) and validation (20%) set.

## 3.4 Model architecture for the reactivity prediction

Various model architectures for reactivity prediction have been
considered, ranging from linear multivariate regression, k-nearest
neighbors (KNN), random forests (RF), graph neural networks
(GNN), XGBoost and FFNN, to geometry-based kernel ridge re-
gression (KRR) models. An in-depth discussion of the character-
istics of each model considered can be found in Section S3 of the
Supporting Information.

## 3.5 Datasets

In total, 3 application domains have been considered to showcase
the ability of our model architectures to efficiently learn barrier
heights in various settings.

In the first application domain, we focused on HAT predic-
tions related to chemical synthesis. We started by analyzing a
dataset consisting of 238 computed reaction profiles for alkoxy
radicals abstracting hydrogens from hydrocarbons and heterosub-
stituted compounds in an acetonitrile solution.⁷⁰ Geometry opti-
mizations and frequency calculations were performed at ωB97X-
D/6-31G(d) level of theory⁸⁴⁻⁸⁷, single-point energy refinements
at UωB97X-D/6-311++G(d,p) level of theory, and solvation was
taken into account through the SMD model.⁸⁸ Since two of the
radicals involved in this dataset contain Cl and F atoms – ele-

ments that are not present in the training data of the surrogate
model – we computed the descriptors for these species explicitly.

Subsequently, we aimed to apply our trained models to an
experimental dataset of rate constants/activation energies for
45 HAT reactions from C(sp³)-H bonds by cumyloxyl radical
(CumO•), compiled by Bietti et al.$^{54}$ Finally, we also considered a
small dataset of experimentally reported selectivities for 6 hydro-
carbons by CH₃O• determined by Zuo et al.$^{50}$ More specifically,
we aimed to reproduce the observed selectivities by inferring the
relative barrier heights associated with individual carbon attack
positions (15 in total).

In a second application, we focused on the prediction of
metabolites of medicinally relevant compounds. As such, we
considered a second dataset consisting of 24 activation en-
ergies for hydrogen atom transfer by the cytochrome P450
enzyme from organic compounds,⁴⁶ computed at B3LYP/6-
311++G(2d,2p)//B3LYP/6-31G(d) (Fe=SVP) level of the-
ory.⁸⁹⁻⁹¹ 5 of the HAT reactions involve compounds with ele-
ments outside of the domain of the training data of our surro-
gate model and hence were discarded in our analysis. The mod-
els trained on the resulting extremely small dataset were sub-
sequently applied to predict activation energies for 6 additional
computed HAT reaction profiles for unseen substrates.

As a third application, we focus on HAT reactions encountered
in atmospheric chemistry. To this end, around 300 HAT reactions
were extracted from RMechDB. RMechDB is a database consisting
of elementary radical reactions recently assembled by Tavakoli et
al..³⁸ In its current form, the dataset is dominated by gas-phase
reactions of isoprene and other organic molecules. No labels, i.e.,
computed or predicted activation energies, are available for the
reactions in this database, and hence we computed reaction pro-
files ourselves with the same high-throughput workflow used to
generate the previously discussed in-house HAT dataset.

# 4 Results and Discussion

## 4.1 Performance of the surrogate model

In Figure 3, the performance of the trained surrogate model on
the held-out test set, i.e., the radicals and closed-shell molecules
present in the BDE-db dataset, is presented. For (relaxed) BDFEs,
atomic charges, spin densities, and buried volumes, excellent cor-
relations between predicted and computed values are obtained
(R² ≥ 0.97); for frozen BDEs, the correlation is slightly worse
(R² = 0.89), but the quantitative trends are still largely retained.

## 4.2 Predicting activation energies for our in-house generated dataset

With a performant surrogate model at our disposal, we turned to
the in-house generated HAT dataset. In Figure 4, the distribution
of the computed activation and reaction energies, as well as the
correlation between both quantities, is presented. Remarkably,
activation and reaction energies only correlate moderately (R² =
0.71), indicating that the Bell-Evans-Polanyi principle⁹²,⁹³ is not
sufficient to fully explain the trends in the activation energies, and
that building a predictive model including multiple descriptors is
warranted for this dataset.
<figure><img src="image_3.png" /><figcaption>Fig. 2 Schematic representation of the surrogate model used for the QM descriptor predictions.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3 Correlation between the computed QM descriptors and those predicted by the surrogate model.</figcaption></figure>

Various baseline models, not making use of the informative representation predicted by our surrogate GNN, were considered as baselines. Due to the relatively small size and inherent diversity present in the constructed HAT dataset, all of these model architectures perform remarkably poorly. RF, KNN and XGBoost models taking structural fingerprints, i.e., differential reaction fingerprints (DRFP), as input, reach a mean absolute error (MAE) of 5.2-5.3 kcal/mol (corresponding to $R^2$ around 0) in 10-fold cross-validation. This is barely any better than the standard deviation (7 kcal/mol), i.e., the hypothetical accuracy that would be obtained for a model that always predicts the mean value for the dataset. The GNN architectures tested – a Weisfeiler-Lehmann (WL) one and Chemprop – perform somewhat better, but the MAE still amounts to ~2.9-3.5 kcal/mol ($R^2$ ~ 0.65). The geometry-based kernel ridge regression model, based on the recently introduced Bond-Based Reaction Representation ($B^2R_i^2$) reaches a similar performance (MAE = 3.55 kcal/mol).⁹⁴

The results quoted above are in great contrast to the results obtained when the predicted informative representation is used as the input for the reactivity model. Now, the MAEs amount to 1.9-2.3 kcal/mol (cf. Table 1) and the corresponding $R^2$ scores amount to 0.80-0.84. There is still some non-linearity present in the relationship between the representation used and the activation energies but only to a limited extent: our best non-linear model architecture, the FFNN, outperforms a multivariate linear regression model by only 0.30 and 0.36 kcal/mol, in terms of MAE and RMSE, respectively (ensembling 4 FFNNs, the error can be reduced further by 0.02 and 0.03 kcal/mol respectively). In line with this observation, we also observe that fairly few data points are needed to reach this accuracy; already with 800 training points, the final accuracy is approached within 0.1 kcal/mol for most of the machine learning model architectures (Figure 5). Reducing the training set size, the multivariate linear model takes over as the most accurate model when fewer than 400 training
<figure><img src="image_5.png" /><figcaption>Fig. 4 (a) Histogram presenting the distribution of the computed activation energies with tunneling corrections included ($\Delta G_{corr}^{\ddagger}$). (b) Histogram representing the distribution of the computed reaction energies ($\Delta G_{rxn}$). (c) Correlation plot between $\Delta G_{corr}^{\ddagger}$ and $\Delta G_{rxn}$.</figcaption></figure>

points are considered.

<table><caption>Table 1 A summary of the performance of some of the model architectures tested on the prediction of $\Delta G^{\ddagger}$ for the HAT dataset.</caption><thead><tr><th>model architecture</th><th>MAE (kcal/mol)</th><th>RMSE (kcal/mol)</th></tr></thead><tbody><tr><td>RF (fingerprint input)</td><td>5.25</td><td>6.96</td></tr><tr><td>XGBoost (fingerprint input)</td><td>5.32</td><td>7.03</td></tr><tr><td>Chemprop (CGR representation)</td><td>2.94</td><td>4.15</td></tr><tr><td>KRR ($B^2R_l^2$ representation)</td><td>3.55</td><td>4.71</td></tr><tr><td>linear model (descriptor input)</td><td>2.28</td><td>3.14</td></tr><tr><td>RF (descriptor input)</td><td>2.10</td><td>3.01</td></tr><tr><td>XGBoost (descriptor input)</td><td>2.24</td><td>3.13</td></tr><tr><td>FFNN (descriptor input)</td><td>1.98</td><td>2.78</td></tr></tbody></table>

In Fig. 6a, a schematic overview of the impact of the individual descriptors of the (predicted) informative representation on the accuracy of the multivariate linear model is provided. Our analysis indicates that each descriptor has at least some predictive power, i.e., they all correlate to some extent with the activation energies, so that the prediction accuracy for the corresponding univariate linear regression models exceeds the standard deviation across the dataset by a non-negligible margin. As could have been expected from the Bell-Evans-Polanyi analysis in Fig. 4, the (relaxed) BDFE values correlate extremely well; a linear regression model based exclusively on this (predicted) descriptor reaches an MAE of 2.84 kcal/mol, which is not dramatically worse than the multivariate model based on all (predicted) descriptors. The descriptor that is most complementary to the BDFE, i.e., the descriptor whose inclusion in the linear model reduces the prediction errors the most, is the atomic charge ($q$), closely followed by the frozen bond dissociation energy ($BDE_{fr}$). On their own, the spin density ($\rho$) values on the radical centers of reactants and products, also have significant predictive power, but they do not improve the model accuracy when the BDFE is already considered as an input descriptor. Note that this is in line with the

point of view discussed above that changes in delocalization are fully encoded in the thermodynamic driving force: the information encoded in the spin densities overlaps with the information gained from consideration of the BDFE values.⁶⁴ Remarkably, and in contrast to the results obtained for the smaller subsets of chemical space considered as part of the applications (vide infra), the inclusion of $V_{bur}$ does not improve the prediction accuracy on this dataset.

## 4.3 Applications

In this final section, we will showcase the usefulness of our approach by considering a set of applications. Particularly, we will exploit the data efficiency of models based on our VB-inspired representation to learn activation energies across some extremely small HAT datasets found in the literature (vide supra).

As indicated above, the first application area focuses on HAT predictions related to chemical synthesis. More specifically, we started by considering a HAT dataset involving alkoxy radicals and hydrocarbons/heterosubstituted compounds.⁷⁰ In previous work, an empirical model based on five computed parameters, respectively describing the electronegativities and delocalization of the reactants and products, as well as the thermodynamic driving force, was proposed.⁵¹ Applying this empirical model to the full dataset, an MAE 0.85 kcal/mol ($R^2=0.84$) is obtained on the training set, and 1.24 kcal/mol ($R^2=0.72$) on the test set.

With the help of our on-the-fly predicted VB-inspired representation, MAEs of 0.99 kcal/mol ($R^2=0.76$) and 1.14 kcal/mol ($R^2=0.72$) are reached, respectively on the training and test set, with a simple multivariate linear model. While not significantly better than the results obtained through the tailored, empirical expression by Liu *et al.*, it is important to note that our predictions can be performed on a millisecond scale, whereas the computed descriptor inputs needed in the empirical model require geometry generation and a succession of tedious DFT calculations.

We also considered the performance of non-linear model ar-
<figure><img src="image_6.png" /><figcaption>Fig. 5 Learning curves for the $\Delta G^{\ddagger}_{corr}$ for several models based in (a) an intermediate VB learned representation and (b) the differential reaction fingerprints. (c) Correlation plot between the computed and predicted activation energies with tunneling corrections $\Delta G^{\ddagger}_{corr}$ for the (in-sample) test sets considered during 10-fold cross-validation for the FFNN with 4 ensembles.</figcaption></figure>

architectures on this small dataset. Training a FFNN from scratch results in an MAE on the test set of 1.45 kcal/mol ($R^2=0.57$). A transfer learning approach was attempted as well, where the parameters from the FFNN model trained on the more extensive and diverse HAT dataset discussed above were used as the initialization of the network for this dataset. Taking this approach, an MAE on the test set of 1.34 kcal/mol ($R^2=0.63$) is reached; still worse than the linear model. Ensembling 4 FFNNs, the model shows a modest improvement: training the ensemble from scratch, an MAE of 1.43 kcal/mol ($R^2=0.58$) is reached, but using a transfer learning approach, an MAE of 1.21 kcal/mol ($R^2=0.69$) is obtained. Remarkably, an RF model based on the predicted representation now actually outperforms all other model architectures, with an MAE on the test set of 1.07 kcal/mol ($R^2=0.80$).

Another remarkable observation is that the predictive power of the respective descriptors in our VB-inspired representation has changed significantly upon switching from the HAT dataset discussed in the previous section to the current one (Figure 6; cf. Section S5.1 in the Supporting Information). As mentioned before, in the diverse dataset constructed as part of this study, the inclusion of the predicted BDFEs reduced the MAE of the linear model the most, and the inclusion of $V_{bur}$ had no observable effect on the model accuracy. In the small alkoxy radical dataset, however, BDFE values are not nearly as predictive. Atomic charges improve the model accuracy the most now, and the inclusion of $V_{bur}$ leads to significant further improvements. These differences can be connected to the more pronounced polar effects in alkoxy radicals, due to the involvement of a more electronegative radical center.⁵⁷,⁹⁵⁻⁹⁷ In more general terms, these findings about descriptor importance underscore that local correlations in small subsections of chemical space are not necessarily transferable and that caution is needed when aiming to draw universal conclusions about reactivity trends based purely on empirical trends.

Next, we turn to the experimental dataset of 45 HAT activation energies, assembled by Bietti and co-workers. Our models, trained on the alkoxy radical dataset, predict the (rescaled) activation energies in this experimental set with acceptable accuracy

($\sim$ 1 kcal/mol; cf. Section S5.2 of the Supporting Information), correctly identifying all of the reaction barriers to be low. However, the coefficients of determination, i.e., the extent to which fine-grained trends across these 45 reactions are recovered, are remarkably poor ($R^2 \sim 0.35$). The main reason for this failure to recover the reactivity trends is the distinct data distribution, i.e., the activation energies in this specific experimental dataset are "outliers" with respect to the reactions present in Ma *et al.'s* dataset, and – to a lesser extent – the fairly low spread in the targets (the standard deviation amounts to only 1.5 kcal/mol). This conjecture can be confirmed by retraining our models on the 45 data points by Bietti *et al.*: in 10-fold cross-validation (i.e., when only 36 in-distribution training points are considered), a slightly better $R^2 = 0.47$ is now reached for our best model, the ensembled FFNN with transfer learning.

Finally, we also attempted to predict the experimentally reported HAT selectivities reported by Zuo *et al.*,⁵⁰ which were subsequently computationally recovered through DFT calculations by Ma and co-workers.⁷⁰ In Figure 7, the performances of the various model architectures are presented. For all models, the predictions are in qualitative agreement with the experimentally observed regioselectivity for at least 4 out of 6 compounds (with our best model – the ensembled FFNN – obtaining a perfect score of 6 out of 6). Also quantitatively, our models reach good accuracy as compared to the computed DFT barriers in this application. The most accurate model is again the FFNN with 4 ensembles (MAE = 0.84 kcal/mol, RMSE = 1.07 kcal/mol, $R^2 = 0.74$), with only one prediction resulting in an error over 2 kcal/mol with respect to the DFT calculations. The RF model, despite being the least successful at qualitatively predicting the most likely HAT position, reaches only a slightly worse quantitative accuracy (MAE = 0.94 kcal/mol, RMSE = 1.37 kcal/mol; $R^2 = 0.64$).

The second application area focuses on P450 metabolism of medicinally relevant compounds. As discussed in the Methodology Section, predictive models were trained on a tiny dataset of 19 HAT activation energies, and subsequently, these models were applied to predict 6 barriers for unseen substrates, in line
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6 Performance on ΔG‡, in terms of mean absolute error (MAE), for multivariate linear models based on a subset of the descriptors. (a) Trained in in-house HAT dataset. (b) Trained in 238 alkoxy radicals abstracting hydrogens.</figcaption></figure>

with the work by Gingrich *et al.*.⁴⁶ We again started with a multi-variate linear regression model, but we considered only (a subset of) the descriptors of the hydrocarbons to reduce collinearity and avoid overfitting (Section S6). With the trained model, based only on the predicted spin density and buried volume, an MAE of 1.12 kcal/mol is obtained for the unseen substrates, in line with the accuracy reached by the predictive model proposed by Gingrich *et al.* (1.15 kcal/mol).⁴⁶ Adding more descriptors, the MAE of the model increases (Section S6).

Subsequently, we also considered our FFNN architecture. Trained from scratch on the same subset of descriptors with two neurons, the model reaches an MAE of 1.04 kcal/mol. With 4 ensembles, the model reaches an MAE of 0.94 kcal/mol.

For the final application, we considered the HAT reactions found in RMechDB (Section S7 of the Supporting Information). In total, we were able to extract 268 reactions, falling within the scope of our surrogate model, from this database; primarily HAT reactions between organic compounds and HO· and HOO· radicals, i.e., the onset of atmospheric oxidation. With our autodE-based workflow for reaction profile computation, we only managed to successfully compute profiles for 73 of these reactions (failure sources are detailed in Section S7).

At first, we tried to directly reproduce the computed activation energies with the predictive models trained on our in-house assembled HAT dataset, but even our best model (ensembled FFNN) reached an inadequate MAE of 4.58 kcal/mol (RMSE = 5.11 kcal/mol, $R^2=0.05$). A reasonable explanation for this poor performance is (once again) the distribution of the data: our in-house generated dataset contained no examples involving HO· or HOO· radicals, and the activation energies for the reactions extracted from RMechDB (mean $ΔG_{corr}^‡$ = 9.10 kcal/mol) tend to be at the lower edge of the activation energy distribution of our training data (cf. Figure 4). Once

more, upon re-training on the 73 computed RMechDB reaction profiles, an outstanding accuracy is recovered; in 10-fold cross-validation, our ensembled FFNN with transfer learning now reaches $R^2$ = 0.88 (MAE = 1.29 kcal/mol, RMSE = 1.51 kcal/mol); our multivariate linear model reaches $R^2$ = 0.88 as well (MAE = 1.20 kcal/mol, RMSE = 1.57 kcal/mol). With both of these fine-tuned models, we made predictions for the remaining reactions extracted from the database (cf. https://github.com/chimie-paristech-CTM/energy predictor HAT/tree/master/reactivity model/RMechDB pred).

The observations for this last application are yet another illustration that setting up a machine learning model for chemical reactivity, able to generalize across application domains, is hard and data-intensive, but that introducing QM-based descriptors in the model architecture – particularly in combination with transfer learning – can remedy this situation to a large extent.

## 5 Conclusions

The poor data-efficiency of generic machine learning models poses a major challenge to their application in chemical reactivity studies, since for most specialized predictive tasks in this domain, only limited training data is available. In this work, we have explored a strategy to make machine learning workflows for reactivity prediction significantly less data-hungry, by making use of an intermediate informative reaction representation. This informative representation was constructed based on a VB analysis for generic HAT reactions. A surrogate model, able to predict (most of) the identified quantum chemical informative descriptors was subsequently built based on a repurposed, publicly available dataset of descriptors and geometries of organic radicals. Combining the surrogate model with a secondary reactivity model, we achieved unequivocally superior performance on a medium-sized, in-house generated, dataset of DFT computed re-
<table><thead><tr><th>Molecule</th><th>Site</th><th>Exp ratio</th><th>DFT</th><th>Linear regression</th><th>Random Forest</th><th>FFNN (1 ensemble)</th><th>FFNN (4 ensembles)</th></tr></thead><tbody><tr><td rowspan="2">1<br>2</td><td>1</td><td>97</td><td>12.47</td><td>12.43</td><td>13.52</td><td>12.94</td><td>12.91</td></tr><tr><td>2</td><td>3</td><td>17.28</td><td>16.57</td><td>17.15</td><td>15.95</td><td>16.81</td></tr><tr><td rowspan="2">1<br>2</td><td>1</td><td>47</td><td>14.34</td><td>13.82</td><td>15.11</td><td>14.46</td><td>14.41</td></tr><tr><td>2</td><td>48</td><td>14.05</td><td>13.69</td><td>14.43</td><td>14.36</td><td>14.07</td></tr><tr><td rowspan="3">1<br>2<br>3</td><td>3</td><td>5</td><td>17.56</td><td>16.29</td><td>17.13</td><td>15.98</td><td>16.55</td></tr><tr><td>1</td><td>-</td><td>17.23</td><td>17.65</td><td>15.84</td><td>15.38</td><td>16.35</td></tr><tr><td>2</td><td>-</td><td>15.62</td><td>18.80</td><td>15.38</td><td>14.46</td><td>15.25</td></tr><tr><td rowspan="3">1<br>2</td><td>3</td><td>1</td><td>13.35</td><td>14.67</td><td>15.56</td><td>14.69</td><td>14.87</td></tr><tr><td>1</td><td>22</td><td>14.97</td><td>13.38</td><td>15.10</td><td>14.66</td><td>14.70</td></tr><tr><td>2</td><td>78</td><td>13.01</td><td>20.30</td><td>15.27</td><td>14.03</td><td>14.58</td></tr><tr><td rowspan="3">1<br>2<br>3</td><td>1</td><td>14</td><td>15.71</td><td>13.97</td><td>15.08</td><td>14.33</td><td>14.60</td></tr><tr><td>2</td><td>77</td><td>12.10</td><td>12.47</td><td>12.97</td><td>13.97</td><td>12.84</td></tr><tr><td>3</td><td>9</td><td>16.87</td><td>16.30</td><td>16.75</td><td>15.92</td><td>16.66</td></tr><tr><td rowspan="2">1<br>2</td><td>1</td><td>1</td><td>10.49</td><td>11.46</td><td>13.36</td><td>13.64</td><td>13.05</td></tr><tr><td>2</td><td>-</td><td>12.84</td><td>12.77</td><td>13.43</td><td>14.39</td><td>14.18</td></tr></tbody></table>

Fig. 7 Prediction of selectivity and ΔG‡ for the compounds experimentally considered by Zuo et al. for the various model architectures trained on the alkoxy radical dataset.

action profiles of diverse HAT reactions, compared to any other machine learning approach we tried. Because of the extreme data efficiency of our descriptor-augmented models, we were able to demonstrate that they can be fine-tuned and applied to small datasets across various reaction conditions, settings and application domains, ranging from regular (liquid phase) synthesis, over metabolism and drug design, to atmospheric chemistry.

The presented modular approach towards reactivity prediction has the potential to significantly expand the scope of machine learning workflows for focused reactivity predictions, as long as the reactivity type is qualitatively understood and can be modeled accurately, and high-quality descriptor datasets are available for surrogate model training. In this context, it is crucial to emphasize the necessity for careful curation of these property datasets, ensuring that the chosen descriptors are pertinent, i.e., relevant, to a diverse array of reaction types.

## 6 Data availability

The main code used to generate the presented results can be found at https://github.com/chimie-paristech-CTM/bde-hat. The code associated with the final reactivity model, i.e., the surrogate model in combination with an FFNN, can be found at https://github.com/chimie-paristech-CTM/energy-predictor_HAT. The (generated/curated) datasets can be downloaded at https://figshare.com/projects/Hydrogen_atom-transfer_reactions/188007 and https://figshare.com/articles/dataset/dataset_surrogate_model_cs/24754341 respectively.

## Author Contributions

JEAR: data curation, software, methodology, formal analysis, visualization, writing. RMN: software, methodology, writing (review), visualization. TS: conceptualization, methodology, formal analysis, writing, supervision, funding acquisition.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

JEAR and TS acknowledge the French National Agency for Research (ANR) for a CPJ grant (ANR-22-CPJ1-0093-01). This work was granted access to the HPC resources of IDRIS under the allocation 2023-100732 granted by GENCI. RMN thanks the Zeno Karl Schindler Foundation for their generous financial support, and Connor Coley and Gisbert Schneider for their guidance during her Master's thesis work.

## Notes and references

1. D. Weininger, *J. Chem. Inf. Comput. Sci.*, 1988, **28**, 31–36.

2. C. W. Coley, R. Barzilay, T. S. Jaakkola, W. H. Green and K. F. Jensen, *ACS Cent. Sci.*, 2017, **3**, 434–443.

3. C. W. Coley, W. Jin, L. Rogers, T. F. Jamison, T. S. Jaakkola, W. H. Green, R. Barzilay and K. F. Jensen, *Chem. Sci.*, 2019, **10**, 370–377.

4. P. Schwaller, T. Laino, T. Gaudin, P. Bolgar, C. A. Hunter, C. Bekas and A. A. Lee, *ACS Cent. Sci.*, 2019, **5**, 1572–1583.
5 M. S. Sigman, K. C. Harper, E. N. Bess and A. Milo, *Acc. Chem. Res.*, 2016, **49**, 1292–1301.

6 A. F. Zahrt, J. J. Henle, B. T. Rose, Y. Wang, W. T. Darrow and S. E. Denmark, *Science*, 2019, **363**, eaau5631.

7 D. T. Ahneman, J. G. Estrada, S. Lin, S. D. Dreher and A. G. Doyle, *Science*, 2018, **360**, 186–190.

8 W. Beker, E. P. Gajewska, T. Badowski and B. A. Grzybowski, *Angew. Chem., Int. Ed.*, 2019, **58**, 4515–4519.

9 K. Jorner, T. Brinck, P.-O. Norrby and D. Buttar, *Chem. Sci.*, 2021, **12**, 1163–1175.

10 X. Li, S.-Q. Zhang, L.-C. Xu and X. Hong, *Angew. Chem., Int. Ed.*, 2020, **59**, 13253–13259.

11 D. F. Nippa, K. Atz, R. Hohler, A. T. Müller, A. Marx, C. Bartelmus, G. Wuitschik, I. Marzuoli, V. Jost, J. Wolfard et al., *Enabling late-stage drug diversification by high-throughput experimentation with geometric deep learning*, 2022.

12 K. Low, M. L. Coote and E. I. Izgorodina, *J. Chem. Theor. Comput.*, 2023, **19**, 1466–1475.

13 K. Low, M. L. Coote and E. I. Izgorodina, *J. Chem. Theor. Comput.*, 2022, **18**, 1607–1618.

14 K. Low, R. Kobayashi and E. I. Izgorodina, *J. Chem. Phys.*, 2020, **153**, 104101.

15 J. G. Estrada, D. T. Ahneman, R. P. Sheridan, S. D. Dreher and A. G. Doyle, *Science*, 2018, **362**, eaat8763.

16 K. Yang, K. Swanson, W. Jin, C. Coley, P. Eiden, H. Gao, A. Guzman-Perez, T. Hopper, B. Kelley, M. Mathea et al., *J. Chem. Inf. Model.*, 2019, **59**, 3370–3388.

17 Z. Liu, L. Lin, Q. Jia, Z. Cheng, Y. Jiang, Y. Guo and J. Ma, *J. Chem. Inf. Model.*, 2021, **61**, 1066–1082.

18 R. Zubatyuk, J. S. Smith, J. Leszczynski and O. Isayev, *Sci. Adv.*, 2019, **5**, eaav6490.

19 K. T. Schütt, H. E. Sauceda, P.-J. Kindermans, A. Tkatchenko and K.-R. Müller, *J. Chem. Phys.*, 2018, **148**, 241722.

20 Y. Guan, C. W. Coley, H. Wu, D. Ranasinghe, E. Heid, T. J. Struble, L. Pattanaik, W. H. Green and K. F. Jensen, *Chem. Sci.*, 2021, **12**, 2198–2208.

21 W. Jin, C. Coley, R. Barzilay and T. Jaakkola, *Advances in neural information processing systems*, 2017, 30.

22 T. Stuyver and C. W. Coley, *J. Chem. Phys.*, 2022, **156**, 084104.

23 A. Gaulton, L. J. Bellis, A. P. Bento, J. Chambers, M. Davies, A. Hersey, Y. Light, S. McGlinchey, D. Michalovich, B. Al-Lazikani et al., *Nucleic Acid Res.*, 2012, **40**, D1100–D1107.

24 *Pistachio* (NextMove Software), https://www.nextmovesoftware.com/pistachio.html.

25 R. Ramakrishnan, P. O. Dral, M. Rupp and O. A. Von Lilienfeld, *Sci. Data*, 2014, **1**, 1–7.

26 M. Nakata and T. Shimazaki, *J. Chem. Inf. Model.*, 2017, **57**, 1300–1308.

27 M. Nakata, T. Shimazaki, M. Hashimoto and T. Maeda, *J. Chem. Inf. Model.*, 2020, **60**, 5891–5899.

28 A. Stuke, C. Kunkel, D. Golze, M. Todorović, J. T. Margraf, K. Reuter, P. Rinke and H. Oberhofer, *Sci. Data*, 2020, **7**, 58.

29 P. C. St. John, Y. Guan, Y. Kim, B. D. Etz, S. Kim and R. S.

Paton, *Sci. Data*, 2020, **7**, 244.

30 C. Isert, K. Atz, J. Jiménez-Luna and G. Schneider, *Sci. Data*, 2022, **9**, 1–11.

31 R. M. Neeser, C. Isert, T. Stuyver, G. Schneider and C. W. Colley, *Chem. Data. Coll.*, 2023, **46**, 101040.

32 K. Schütt, P.-J. Kindermans, H. E. Sauceda Felix, S. Chmiela, A. Tkatchenko and K.-R. Müller, *Advances in neural information processing systems*, 2017, 30.

33 J. Gasteiger, S. Giri, J. T. Margraf and S. Günnemann, *Fast and uncertainty-aware directional message passing for non-equilibrium molecules*, 2020.

34 G. A. Pinheiro, J. Mucelini, M. D. Soares, R. C. Prati, J. L. Da Silva and M. G. Quiles, *J. Phys. Chem. A*, 2020, **124**, 9854–9866.

35 N. De Cao and T. Kipf, *MolGAN: An implicit generative model for small molecular graphs*, 2018.

36 Y. Bian and X.-Q. Xie, *J. Mol. Model.*, 2021, **27**, 1–18.

37 J. Westermayr, J. Gilkes, R. Barrett and R. J. Maurer, *Nat. Comput. Sci.*, 2023, **3**, 139–148.

38 M. Tavakoli, Y. T. T. Chiu, P. Baldi, A. M. Carlton and D. Van Vranken, *J. Chem. Inf. Model.*, 2023, **63**, 1114–1123.

39 Y. Zhang, J. Yu, H. Song and M. Yang, *Structure-Based Reaction Descriptors for Predicting Rate Constants by Machine Learning: Application to Hydrogen Abstraction from Alkanes by CH3/H/O Radicals*, 2023.

40 J. Yu, D. Shan, H. Song and M. Yang, *Fuel*, 2022, **322**, 124150.

41 W. Nie, C. Yan, L. Yang, P. Roldin, Y. Liu, A. L. Vogel, U. Molteni, D. Stolzenburg, H. Finkenzeller, A. Amorim et al., *Nat. Commun.*, 2023, **14**, 3347.

42 L. De Keer, P. Van Steenberge, M.-F. Reyniers, G. Gryn'Ova, H. M. Aitken and M. L. Coote, *Polym. Chem.*, 2022, **13**, 3304–3314.

43 G. Gryn'ova, J. L. Hodgson and M. L. Coote, *Org. Biomol. Chem.*, 2011, **9**, 480–490.

44 A. Grinberg Dana, H. Wu, D. S. Ranasinghe, F. C. Pickard IV, G. P. Wood, T. Zelesky, G. W. Sluggett, J. Mustakis and W. H. Green, *Mol. Pharmaceutics*, 2021, **18**, 3037–3049.

45 H. Wu, A. Grinberg Dana, D. S. Ranasinghe, F. C. Pickard IV, G. P. Wood, T. Zelesky, G. W. Sluggett, J. Mustakis and W. H. Green, *Mol. Pharmaceutics*, 2022, **19**, 1526–1539.

46 P. W. Gingrich, J. B. Siegel and D. J. Tantillo, *Chemistry Methods*, 2022, **2**, e202100108.

47 H. Matsubara, T. Kawamoto, T. Fukuyama and I. Ryu, *Acc. Chem. Res.*, 2018, **51**, 2023–2035.

48 S. Sarkar, K. P. S. Cheung and V. Gevorgyan, *Chem. Sci.*, 2020, **11**, 12974–12993.

49 L. Capaldo and D. Ravelli, *Eur. J. Org. Chem.*, 2017, **2017**, 2056–2071.

50 Q. An, Z. Wang, Y. Chen, X. Wang, K. Zhang, H. Pan, W. Liu and Z. Zuo, *J. Am. Chem. Soc.*, 2020, **142**, 6216–6226.

51 F. Liu, S. Ma, Z. Lu, A. Nangia, M. Duan, Y. Yu, G. Xu, Y. Mei, M. Bietti and K. Houk, *J. Am. Chem. Soc.*, 2022, **144**, 6802–6812.

52 L.-C. Yang, X. Li, S.-Q. Zhang and X. Hong, *Org. Chem. Front.*,
2021, **8**, 6187–6195.

53 Y. Sun, F. Liu, J. N. Sanders and K. Houk, *J. Org. Chem.*, 2023, **88**, 12668–12676.

54 M. Salamone, M. Galeotti, E. Romero-Montalvo, J. A. van den B. D. Groff, J. M. Mayer, G. A. DiLabio and M. Bietti, *J. Am. Chem. Soc.*, 2021, **143**, 11759–11776.

55 M. Finn, R. Friedline, N. K. Suleman, C. J. Wohl and J. M. Tanko, *J. Am. Chem. Soc.*, 2004, **126**, 7578–7584.

56 K. R. Korzekwa, J. P. Jones and J. R. Gillette, *J. Am. Chem. Soc.*, 1990, **112**, 7042–7046.

57 B. D. Groff, B. Koronkiewicz and J. M. Mayer, *The Journal of Organic Chemistry*, 2023.

58 S. Shaik and A. Shurki, *Angew. Chem., Int. Ed.*, 1999, **38**, 586–625.

59 S. S. Shaik and P. C. Hiberty, *A chemist's guide to valence bond theory*, John Wiley & Sons, 2007.

60 S. Shaik, D. Kumar and S. P. De Visser, *J. Am. Chem. Soc.*, 2008, **130**, 10128–10140.

61 W. Lai, C. Li, H. Chen and S. Shaik, *Angew. Chem., Int. Ed.*, 2012, **51**, 5556–5578.

62 L. Song, W. Wu, P. C. Hiberty and S. Shaik, *Chem. Eur. J.*, 2006, **12**, 7458–7466.

63 G. N. Sastry and S. Shaik, *J. Am. Chem. Soc.*, 1998, **120**, 2131–2145.

64 T. Stuyver, F. De Proft, P. Geerlings and S. Shaik, *J. Am. Chem. Soc.*, 2020, **142**, 10102–10113.

65 T. Stuyver, D. Danovich, F. De Proft and S. Shaik, *J. Am. Chem. Soc.*, 2019, **141**, 9719–9730.

66 J. Joy, T. Stuyver and S. Shaik, *J. Am. Chem. Soc.*, 2020, **142**, 3836–3850.

67 J. Li, S. Zhou, J. Zhang, M. Schlangen, T. Weiske, D. Usharani, S. Shaik and H. Schwarz, *J. Am. Chem. Soc.*, 2016, **138**, 7973–7981.

68 T. Stuyver and S. Shaik, *J. Am. Chem. Soc.*, 2020, **142**, 20002–20013.

69 L. C. Gallegos, G. Luchini, P. C. St. John, S. Kim and R. S. Paton, *Acc. Chem. Res.*, 2021, **54**, 827–836.

70 S. Ma, S. Wang, J. Cao and F. Liu, *ACS Oomega*, 2022, **7**, 34858–34867.

71 K. Jorner, *morfeus: Molecular features for machine learning*, 2022, https://github.com/digital-chemistry-laboratory/morfeus.

72 S. S. SV, P. C. S. John and R. S. Paton, *Chem. Sci.*, 2021, **12**, 13158–13166.

73 S. S. S. Vejaykummar, Y. Kim, S. Kim, P. S. John and R. Paton, *Expansion of Bond Dissociation Prediction with Machine Learning to Medicinally and Environmentally Relevant Chemical Space*, 2023.

74 T. A. Young, J. J. Silcock, A. J. Sterling and F. Duarte, *Angew. Chem., Int. Ed.*, 2021, **133**, 4312–4320.

75 M. Frisch, G. Trucks, H. Schlegel, G. Scuseria, M. Robb, J. Cheeseman, G. Scalmani, V. Barone, G. Petersson, H. Nakatsuji et al., *Gaussian 16*, 2016.

76 T. Stuyver, K. Jorner and C. W. Coley, *Sci. Data*, 2023, **10**, 66.

77 D. Probst, P. Schwaller and J.-L. Reymond, *Digit. Discov.*, 2022, **1**, 91–97.

78 N. Casetti, J. E. Alfonso-Ramos, C. W. Coley and T. Stuyver, *Chem. Eur. J.*, 2023, e202301957.

79 Y. Zhao and D. G. Truhlar, *Theor. Chem. Acc.*, 2008, **120**, 215–241.

80 A. Schäfer, H. Horn and R. Ahlrichs, *J. Chem. Phys.*, 1992, **97**, 2571–2577.

81 A. Schäfer, C. Huber and R. Ahlrichs, *J. Chem. Phys.*, 1994, **100**, 5829–5835.

82 P. C. St. John, Y. Guan, Y. Kim, S. Kim and R. S. Paton, *Nat. Commun.*, 2020, **11**, 2328.

83 B. C. Garrett and D. G. Truhlar, *J. Phys. Chem.*, 1979, **83**, 2921–2926.

84 J.-D. Chai and M. Head-Gordon, *Phys. Chem. Chem. Phys.*, 2008, **10**, 6615–6620.

85 J.-D. Chai and M. Head-Gordon, *J. Chem. Phys.*, 2008, **128**, 084106.

86 P. C. Hariharan and J. A. Pople, *Theor. Chim. Acta*, 1973, **28**, 213–222.

87 M. M. Francl, W. J. Pietro, W. J. Hehre, J. S. Binkley, M. S. Gordon, D. J. DeFrees and J. A. Pople, *J. Chem. Phys.*, 1982, **77**, 3654–3665.

88 A. V. Marenich, C. J. Cramer and D. G. Truhlar, *J. Phys. Chem. B*, 2009, **113**, 6378–6396.

89 A. D. Becke, *Phys. Rev. A*, 1988, **38**, 3098.

90 C. Lee, W. Yang and R. G. Parr, *Phys. Rev. B*, 1988, **37**, 785.

91 P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, *J. Phys. Chem.*, 1994, **98**, 11623–11627.

92 R. P. Bell, *Proc. R. Soc. A: Math. Phys. Eng. Sci.*, 1936, **154**, 414–429.

93 M. Evans and M. Polanyi, *Trans. Farad. Soc.*, 1936, **32**, 1333–1360.

94 P. van Gerwen, A. Fabrizio, M. D. Wodrich and C. Corminboeuf, *Mach. Learn.: Sci. Technol.*, 2022, **3**, 045005.

95 M. Salamone, F. Basili and M. Bietti, *The Journal of Organic Chemistry*, 2015, **80**, 3643–3650.

96 J. M. Tedder, *Angewandte Chemie International Edition in English*, 1982, **21**, 401–410.

97 A. Ruffoni, R. C. Mykura, M. Bietti and D. Leonori, *Nature Synthesis*, 2022, **1**, 682–695.