3DReact: Geometric deep learning
for chemical reactions

Abstract

Geometric deep learning models, which incorporate the relevant molecular symmetries within the neural network architecture, have considerably improved the accuracy and data efficiency of predictions of molecular properties. Building on this success, we introduce 3DReact, a geometric deep learning model to predict reaction properties from three-dimensional structures of reactants and products. We demonstrate that the invariant version of the model is sufficient for existing reaction datasets. We illustrate its competitive performance on the prediction of activation barriers on the GDB7-22-TS, Cyclo-23-TS and Proparg-21-TS datasets in different atom-mapping regimes. We show that, compared to existing models for reaction property prediction, 3DReact offers a flexible framework that exploits atom-mapping information, if available, as well as geometries of reactants and products (in an invariant or equivariant fashion). Accordingly, it performs systematically well across different datasets, atom-mapping regimes, as well as both interpolation and extrapolation tasks.

These authors contributed equally to this work. National Center for Competence in Research – Catalysis (NCCR-Catalysis), École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland These authors contributed equally to this work. Learning & Adaptive Systems Group, Department of Computer Science, ETH Zurich, 8092 Zurich, Switzerland Learning & Adaptive Systems Group, Department of Computer Science, ETH Zurich, 8092 Zurich, Switzerland National Center for Competence in Research – Catalysis (NCCR-Catalysis), École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland Learning & Adaptive Systems Group, Department of Computer Science, ETH Zurich, 8092 Zurich, Switzerland National Center for Competence in Research – Catalysis (NCCR-Catalysis), École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland

1 Introduction

Physics-inspired representations that take as input the three-dimensional structure^{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13} (as well as, in some cases, electronic structure^{14, 15, 16, 17}) of molecules and transform it into a fixed-length vector, while respecting known physical laws, have a rich history in molecular property prediction.^{2, 8, 4, 6, 7, 13, 18, 19, 20, 5, 1, 21, 3, 9, 10, 22, 23, 24, 25, 12, 26, 27, 28, 29, 30} Common desiderata^{31, 32, 33, 34} for high-performing representations are (i) smoothness, (ii) encoding of the appropriate symmetries to permutations, rotations and translations,^{24, 35} (iii) completeness and (iv) additivity to allow for extrapolation to larger systems. Such fingerprints,^{2, 4, 3, 24, 6, 7, 8, 5, 11, 12, 13} being rooted in fundamental principles, are designed to be property-independent: a single representation can be constructed for a molecule to predict any quantum-chemical target. This is analogous to the molecular Hamiltonian, which specifies the energy and all other properties of a system as a function of atoms’ types and positions in three-dimensional space (assuming the molecules are charge neutral and singlets). These representations are typically used in combination with kernel models due to their data efficiency, ability to deal with high-dimensional feature vectors, and interpretability of the similarity kernel.^{2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 13, 31, 32, 33} Early works showed that combining such representations^{2, 4, 6, 8, 36} with simple feed-forward neural networks instead of kernel models did not necessarily led to better performance.^{37, 38}

More recently, end-to-end neural networks have been proposed that learn the representation as part of the (supervised) training process,^{39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61} based on similar principles to the aforementioned physics-inspired representations: they take as input a three-dimensional structure, as well as in some cases charge and spin information.^{46, 51, 52, 53} The network may be invariant or equivariant to rotations and translations of the input molecules. The former is typically achieved by operating on distances between atoms,^{39, 40, 42} and the latter by operating on relative position vectors and angular information processed by rotationally-equivariant convolutional layers.^{46, 49, 50, 62, 59, 45, 43, 41, 44, 48, 54, 55, 56, 57, 58} Equivariant models are naturally suited to predict vectorial^{44, 49, 62, 43, 59, 48, 45} or higher order tensorial^{59, 52, 54, 55, 63} properties. They have also been demonstrated to exhibit improved data efficiency and generalization capabilities compared to their invariant counterparts on predictions of scalar properties,⁴³ albeit at a higher computational cost. Nevertheless, given an expressive enough architecture (i.e., using higher-order messages^{41, 62, 64, 65, 56, 66, 67} and/or enough convolutional layers^{43, 52, 56}), invariant models are sufficient for many property prediction tasks.⁵⁶

Despite these advances for molecular property prediction, the prediction of computed reaction properties (principally, reaction barriers^{68, 36, 38, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82}) is still in its infancy.⁸³ Machine learning approaches span from utilizing simple two-dimensional fingerprints of reaction components^{84, 85} (reactants and products) to physical-organic descriptors^{86, 87, 76, 88, 89, 75, 90, 91, 92, 93, 94, 95, 80, 96, 97, 98, 82}, or electronic structure-inspired features⁹⁹, to transformer models^{100, 101} adapted for regression,¹⁰² and 2D graph-based approaches^{71, 103, 70, 81}. The latter, particularly the ChemProp model,^{71, 103} are often best-in-class in predicting reaction properties.¹⁰³ It has been shown³⁸ that these models achieve their impressive performance by exploiting atom-mapping information,^{104, 105, 106, 107} which provide information analogous to the reaction mechanism.

Another category of reaction fingerprints arises from discretization of physically-inspired functions^{2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 13} constructed using a cheap estimate of the transition state (TS) structure⁷³ or rather the structures of the reaction components^{36, 69, 74} The SLATM_d representation^{69, 36} in particular has been shown³⁸ to yield accurate predictions of reaction barriers, particularly for datasets^{108, 69} relying on subtle changes in the geometry of reactants and/or products. End-to-end models based on three-dimensional structures of reactants and products have also recently emerged.^{99, 72, 109} In a different vein, several works^{110, 111, 112, 113, 114, 115} aim to directly predict the TS structure, which together with the reactant structure gives the reaction barrier. These approaches lie outside the scope of the property-prediction focus here.

Due to the diversity of challenges posed by different reaction datasets, neither atom-mapping-based models nor 3D-geometry-based models achieve consistently better performance on reaction property prediction tasks.³⁸ To date, no model has been proposed that can incorporate both chemical (atom-maps) and physical (geometry) priors. To address this gap, we introduce 3DReact, a geometric deep learning model that encodes both the three-dimensional structures of reactants and products as well as atom-mapping information or proxies thereof to predict properties of chemical reactions (showcased here for activation energies).

We demonstrate the performance of 3DReact on three datasets of reaction barriers: GDB7-22-TS,¹¹⁶ Cyclo-23-TS,¹¹⁷ and Proparg-21-TS.^{108, 69} As discussed in previous works,³⁸ these datasets present a myriad of challenges for ML models, from the dependence on chemical information¹¹⁶ to the distinction of subtle changes in configurations.^{108, 69} We show that, compared to state-of-the-art models for reaction property prediction,^{103, 36} 3DReact offers accurate and reliable performance across different datasets as well as atom-mapping regimes, reduced dependence on the quality of three-dimensional geometries, and stable extrapolation behavior.

2 Architecture

Figure 1: Architecture of 3DReact. Molecules pass through independent symmetry-adapted (invariant or equivariant) channels (green and orange). These are combined to yield a reaction representation (blue) which is used to predict a reaction property, such as the activation energy (red dot).

3DReact is built from $\mathrm{O}(3)$-equivariant convolutional networks over point clouds as implemented in e3nn.¹¹⁸ Specifically, we use the tensor field network architecture⁴⁷ for molecular components as in Corso et al. ¹¹⁹ While the architecture is equivariant by default, it can easily be made invariant (vide infra). The geometries of molecules constituting reactants and products of each reaction are passed through separate channels, detailed in Section 2.1. They are then combined to eventually predict a reaction property, such as the activation energy, as detailed in Section 2.2. The overall architecture is summarized in Figure 1.

2.1 Symmetry-adapted molecular channels

A molecule with $N_{\mathrm{at}}$ atoms is represented as a distance-based graph where nodes describe atoms and edges describe bonds. Instead of explicitly using connectivity information, the “bonds” of atom $a$ are formed with all the neighboring $\operatorname{Neigh}(a)$ atoms within the cutoff $r_{\max}$. Initial scalar bond (edge) features $\{\mathrm{\mathbf{e}}^{(0)}_{ab}\}$ between atoms $a$ and $b$, as well as spherical harmonics filters $\{\mathrm{\mathbf{z}}_{ab}\}$, are computed from internal coordinates, as detailed in equations LABEL:S-eq:initial-edge-features-start–LABEL:S-eq:initial-edge-features-end. The atom (node) features $\{\mathrm{\mathbf{x}}^{(0)}_{a}\}$ are initialized with $n_{f}{=}16$ cheminformatics descriptors computed with RDKit.¹²⁰ These include atomic number, chirality tag (unspecified, tetrahedral, or other, including octahedral, square planar, allene-type), number of directly-bonded neighbors, number of rings, implicit valence, formal charge, number of attached hydrogens, number of unpaired electrons, hybridization, aromaticity, and presence in rings of specified sizes from $3$ to $7$. This choice is inspired by EquiBind¹²¹ and DiffDock¹¹⁹ and is in line with the improved ⁷² features used for 2D-based methods.

The initial node and edge features pass through embeddings to give $\{\mathrm{\mathbf{x}}^{(1)}_{a}\}$ and $\{\mathrm{\mathbf{e}}_{ab}\}$ respectively, the former are then updated by $n_{\mathrm{conv}}\in\{2,3\}$ equivariant convolutional layers. Each layer is a fully-connected weighted tensor product, as defined in e3nn¹¹⁸. Equations LABEL:S-eq:convs-start–LABEL:S-eq:convs-end describe the equivariant operations performed by the network (see Section LABEL:S-sec:molecular_channels for mathematical details). The network with equivariant molecular components as described is referred to as EquiReact, where its invariant counterpart InReact uses only the $\ell=0$ (scalar) spherical harmonics to construct the convolution filters. The output of the molecular channels is the local molecular representation $\mathrm{\mathbf{X}}\in\mathbb{R}^{N_{\mathrm{at}}\times D}$ corresponding to $N_{\mathrm{at}}$ atoms associated with $D$ features. Depending on the sum_mode hyperparameter, it is constructed either from the node features (node mode) or both node and edge features (both mode).

Inspired by the ChemProp model,^{71, 103} we added an option to exclude hydrogen atoms as nodes when constructing the graph. The only information about hydrogens is then contained in the initial edge features of heavy atoms.

2.2 Combining molecules for reactions

Once atom-wise molecular representations $\mathrm{\mathbf{X}}$ are learned for reactant and product molecules, they must be combined to form a reaction representation $\mathrm{\mathbf{X}}_{\mathrm{rxn}}$.

Figure 2: Scheme illustrating how the reactant (green) and product (orange) representations are combined to form a reaction representation (blue) and eventually predict the target property (red dot) using a multilayer perceptron (mlp). $\sum$ refers to the summation over atom-wise environments. Oblong rectangles and squares represent vectors and scalars, respectively.

For certain datasets, atom-mapping information is available, which correlates individual atoms in reactant molecules to individual atoms in product molecules according to the reaction mechanism. In this setting, the representations $\mathrm{\mathbf{X}}_{\mathrm{reactant}}$ and $\mathrm{\mathbf{X}}_{\mathrm{product}}$ are re-ordered such that the local representation vectors correspond to the same atom in reactants and products. Depending on the combine_mode hyperparameter, either a difference is taken between products’ and reactants’ atom representations, or they are summed, averaged, or passed through a multilayer perceptron (MLP). Thus, the local reaction representation $\mathrm{\mathbf{X}}_{\mathrm{rxn}}$ consists of vectors reflecting how the environment changes in the reaction for each atom. We will refer to this variant of the model, which uses atom-mapping information, as 3DReact_M. While the current model is unable to treat unbalanced reactions (where there are additional atoms on the left- or right-hand side of the reaction equation), its modification in the spirit of ChemProp^{71, 103} is straightforward.

With the reaction representation at hand, predictions are made in the so-called vector or energy modes. In vector mode, the atomic vectors constituting the reaction representation $\mathrm{\mathbf{X}}_{\mathrm{rxn}}$ are initially passed through an MLP to introduce nonlinearity and then summed up to form a global reaction representation vector $\mathrm{\mathbf{\bar{X}}}_{\mathrm{rxn}}$. The target is then learned using an MLP. This model pipeline is illustrated in Figure 2a. In energy mode, on the other hand, the local reaction representations are used to learn atomic contributions to the target (Figure 2b). While performing worse in general, in some cases this mode yields the best predictions (see Section 3.2.1).

Atom-mapping provides static information, analogous to a reaction mechanism, to link atoms in reactants to atoms in products. While highly informative, and thought to be critical to the performance of 2D-graph-based models, ^{71, 103, 72, 70, 81} accurate atom-maps are not available for all reaction datasets.^{38, 104, 105} To circumvent the need for atom-mapping, but mimic its role in exchanging information between reactants and products, other approaches dynamically (i.e., in a learnable fashion) exchange information between molecular representations. For example, RXNMapper¹⁰⁷ is a neural network that learns atom-mappings within the larger self-supervised task of predicting the randomly masked parts in a reaction sequence, using one head of a multi-head transformer architecture. EquiBind,¹²¹ a neural network that predicts the rotation and translation of a ligand to a protein, contains a cross-attention module between ligand and receptor. The latter inspires our surrogate for atom-mapping: 3DReact_X also uses cross-attention between reactants and products to link their atom indices (Section LABEL:S-sec:cross). The re-ordered representations of reactants and products are combined as for the case of atom-mapped reactions (Figures 2a and 2b). We note that other algorithms could also have been used to exchange information between reactants and products, for example in the form of message passing or equivariant attention.^{57, 122}

3DReact also has a simple “no mapping” variant, called 3DReact_S, which does not rely on atom-mapping, nor a surrogate cross-attention module. In vector mode (Figure 2c), the atomic components of molecular representations $\mathrm{\mathbf{X}}_{\mathrm{reactant}}$ and $\mathrm{\mathbf{X}}_{\mathrm{product}}$ are summed up to obtain global vectors $\mathrm{\mathbf{\bar{X}}}_{\mathrm{reactant}}$ and $\mathrm{\mathbf{\bar{X}}}_{\mathrm{product}}$, respectively. Then they are combined, according to the combine_mode parameter, to form a reaction vector $\mathrm{\mathbf{\bar{X}}}_{\mathrm{rxn}}$ which is used to learn the target with an MLP. In energy mode (Figure 2d) individual atomic representations are used to learn their contributions to the quasi-molecular energies of reactants and products, which are later combined (according to the combine_mode parameter) to predict the target. In most cases, this simpler model out-performs 3DReact_X (vide infra).

3 Results and Discussion

The performance of 3DReact is reported for three diverse datasets (the GDB7-22-TS,¹¹⁶ Cyclo-23-TS ¹¹⁷ and Proparg-21-TS ^{108, 69}) using both random and extrapolative splits. For details on the datasets, refer to Section 5.1. For details on the extrapolation splits, see Section 5.2.

Models are run in three atom-mapping regimes: (i) with high-quality maps (“True”) derived from the TS structures or heuristic rules;^{117, 123, 116, 71, 106} (ii) with atom-maps obtained using the open-source RXNMapper¹⁰⁷ (“RXNMapper”); and (iii) without any atom-mapping information at all (“None”). As discussed in recent work,^{124, 38} previously developed graph-based models for reaction property prediction^{71, 72, 70, 96, 97} including ChemProp^{71, 103} reported prediction errors only in the “True” atom-mapping regime. The “RXNMapper” regime is important for cases where the reaction mechanism is not known and atom-mapping using heuristic rules is impossible. The “None” regime is critical for all chemistry that falls outside the realm of organic chemistry captured in the patents¹²⁵ that RXNMapper¹⁰⁷ is trained on.

The atom-mapping-based model 3DReact_M is used in the “True” and “RXNMapper” regimes. In the “None” regime, 3DReact_X and 3DReact_S were tested. 3DReact_S consistently outperformed 3DReact_X, so we include only 3DReact_S and refer the reader to Section LABEL:S-sec:cross for their comparison.

3.1 Equivariance vs. invariance

Table 1 compares the relative performance of the invariant (InReact) and the equivariant (EquiReact) implementations of 3DReact with the learning curves of the two models presented in Figure 3. Previous studies^{43, 56} demonstrated that the equivariant models showed superior extrapolation capabilities on predictions of energies and forces, as well as steeper and shifted learning curves in force prediction tasks. Instead, we find that InReact and EquiReact are practically indistinguishable for the present chemical reaction tasks.

Dataset (property, units)	Atom-mapping regime	InReact	EquiReact
Random splits
GDB7-22-TS ($\Delta E^{\ddagger}$, kcal/mol)	True	$4.93\pm 0.18$	$4.93\pm 0.15$
	RXNMapper	$6.03\pm 0.26$	$6.05\pm 0.25$
	None	$6.56\pm 0.26$	$6.53\pm 0.28$
Cyclo-23-TS ($\Delta G^{\ddagger}$, kcal/mol)	True	$2.39\pm 0.08$	$2.30\pm 0.09$
	RXNMapper	$2.37\pm 0.07$	$2.35\pm 0.12$
	None	$2.39\pm 0.05$	$2.31\pm 0.09$
Proparg-21-TS ($\Delta E^{\ddagger}$, kcal/mol)	True	$0.33\pm 0.07$	$0.31\pm 0.05$
	None	$0.34\pm 0.06$	$0.31\pm 0.06$
Scaffold splits
GDB7-22-TS ($\Delta E^{\ddagger}$, kcal/mol)	True	$7.8\pm 0.7$	$7.8\pm 0.8$
	RXNMapper	$9.2\pm 0.8$	$9.1\pm 0.8$
	None	$10.1\pm 0.9$	$10.0\pm 0.9$
Cyclo-23-TS ($\Delta G^{\ddagger}$, kcal/mol)	True	$2.79\pm 0.18$	$2.72\pm 0.18$
	RXNMapper	$2.77\pm 0.22$	$2.71\pm 0.23$
	None	$2.76\pm 0.22$	$2.72\pm 0.19$
Proparg-21-TS ($\Delta E^{\ddagger}$, kcal/mol)	True	$0.44\pm 0.11$	$0.40\pm 0.08$
	None	$0.45\pm 0.10$	$0.41\pm 0.09$

Table 1: Performance as measured in mean absolute errors (MAEs) of predictions of 3DReact (InReact vs. EquiReact). 3DReact_M is used for the “True” and “RXNMapper” regimes, and 3DReact_S is used for the “None” regime. MAEs are averaged over 10 folds of 80/10/10 splits (training/validation/test) and reported together with standard deviations across folds.

Figure 3: Learning curves for InReact and EquiReact in the “True” atom-mapping regime. Each point shows mean absolute error (MAE), averaged over 10 folds of 80/10/10 splits (for training set fraction $<0.8$, the corresponding subset of the “full” training set is used), and error bars indicate standard deviations across folds.

Figure 4: Top: Reactant and product of a toy reaction: two homometric structures (a) and (b) with atom labels, atom coordinates (Å), and interatomic distances (Å). “Bonds” of the same length are of the same color. Bottom: output after 5 epochs of the invariant (c,d,e) / equivariant (f,g) molecular channels for each atom with different radial cutoffs $r_{\max}$ and number of convolutional layers $n_{\mathrm{conv}}$. Within each subfigure, atomic representations indistinguishable up to shown digits are marked by the same color.

We find that the datasets studied herein do not benefit from the inclusion of equivariant features for molecules. Yet, Figure 4 illustrates that a hypothetical reaction involving conversion between homometric structures of He4,¹²⁶ which is mostly characterized by angle changes, clearly benefits from equivariant molecular features. In the reactant (Figure 4a), all atoms are identical and lead to the same learned representation. In the product (Figure 4b), only atoms B2 and B3 have identical environments, different from A1–4. Atoms B2–3 have the same distances $r$ to the three neighbors, as in A1–4. Thus, InReact, which uses only interatomic distances, yields very close representations for these atoms (Figure 4c). Still, in each convolutional layer, atoms B2 and B3 receive information from B1 and B4, and with increase of $n_{\mathrm{conv}}$ the difference in the representations of B2–3 and A1–A4 becomes more apparent (Figure 4d). However, with smaller radial cutoff $r_{\max}=1.9$, atoms B1–B3 and A1–4 become indistinguishable for any number of layers (Figure 4e). On the other hand, EquiReact, which uses explicit angular information from the spherical harmonics filters, clearly distinguishes all non-equivalent atoms in both cases already for $n_{\mathrm{conv}}=2$ (Figure 4f,g).

While this is a toy example, it illustrates that transformations consisting of changes in angles rather than in bond lengths are better described using EquiReact. In general, the currently available reaction datasets do not pose sufficient challenge to allow distinguishing InReact and EquiReact. For the datasets studied in this work, InReact is sufficient and is the model variation used throughout as 3DReact.

3.2 Benchmark studies

3DReact is compared to previously best baseline models:³⁸ ChemProp,^{71, 103} a graph neural network that uses atom-mapped SMILES to construct a CGR, and the 3D-structure-based SLATM⁸ fingerprint adapted to reactions by taking the difference between product and reactant fingerprints (SLATM_d),³⁶ combined with KRR models (SLATM_d+KRR).

Note that both 3DReact and ChemProp are run without explicit H atoms, for two reasons. First, hydrogen atoms are not always mapped in the “True” and “RXNMapper” regimes, since they are usually implicit in SMILES strings. Second, there is no consistent improvement in including H atoms in the models (Table LABEL:S-tab:models-with-withoutH). SLATM_d, built directly from molecular coordinates without using SMILES strings, does however incorporate H atoms by default. For further discussion refer to Section LABEL:S-sec:hydrogens.

3.2.1 Random splits

Performance as measured in mean absolute errors (MAEs) is illustrated in Table 2 for random splits of each dataset, demonstrating the models’ interpolative capabilities. For the equivalent results with root mean squared errors (RMSEs), consult Section LABEL:S-sec:rmse.

Dataset (property, units)	Atom-mapping regime	ChemProp	SLATMd+KRR	3DReact
GDB7-22-TS ($\Delta E^{\ddagger}$, kcal/mol)	True	$\bf 4.35\pm 0.15$	—	$4.93\pm 0.18$
	RXNMapper	$\bf 5.69\pm 0.17$	—	$\bf 6.03\pm 0.26$
	None	$9.04\pm 0.21$	$\bf 6.89\pm 0.20$	$\bf 6.56\pm 0.26$
Cyclo-23-TS ($\Delta G^{\ddagger}$, kcal/mol)	True	$2.69\pm 0.10$	—	$\bf 2.39\pm 0.08$
	RXNMapper	$2.71\pm 0.07$	—	$\bf 2.37\pm 0.07$
	None	$2.71\pm 0.12$	$2.65\pm 0.08$	$\bf 2.39\pm 0.05$
Proparg-21-TS ($\Delta E^{\ddagger}$, kcal/mol)	True	$1.53\pm 0.14$	—	$\bf 0.33\pm 0.07$
	None	$1.56\pm 0.16$	$\bf 0.33\pm 0.04$	$\bf 0.34\pm 0.06$

Table 2: Performance as measured in mean absolute errors (MAEs) of predictions of 3DReact vs. state-of-the-art baselines ChemProp and SLATM_d+KRR. All datasets are compared in three atom-mapping regimes: “True”, “RXNMapper” and “None”, except for the Proparg-21-TS set, where RXNMapper cannot map the reaction SMILES. MAEs are averaged over 10 folds of random 80/10/10 splits (training/validation/test) and reported together with standard deviations across folds. The lowest errors for each regime and dataset are highlighted in bold, if statistically relevant.

The GDB7-22-TS dataset is distinct from the other two in that it includes variations in the reaction class (and mechanism), thereby showing a greater dependence on the existence and quality of atom-mapping information in the models. It has already been observed³⁸ for ChemProp that there is stark hierarchy in the predictions from the “True” to “RXNMapper” to “None” regimes.

In the “True” regime, 3DReact does not improve predictive capabilities over the ChemProp model for the GDB7-22-TS set. This points to the importance of the chemical diversity in this dataset, where knowledge of the reaction mechanism (in the form of atom-maps) is sufficient information to predict the reaction barriers without information about the geometries of reactants and products. However, as previously discussed,³⁸ “True” maps are an unrealistic scenario for most datasets. Moving to the “RXNMapper” regime, 3DReact and ChemProp already agree within standard deviations. This highlights that for practical-quality maps, 3DReact is amongst the best models for this dataset. In the “None” regime, 3DReact outperforms ChemProp by more than $2\text{\,}\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{/}\mathrm{m}\mathrm{o}\mathrm{l}$.

SLATM_d+KRR results in similar performance to 3DReact for the GDB7-22-TS set. The SLATM_d representation also constructs features from 3D coordinates of the reactants and products using invariant functions, and is therefore more fundamentally similar to 3DReact than ChemProp. Nevertheless, since 3DReact allows for the inclusion of atom-mapping information, predictions are improved in the mapped regimes compared to SLATM_d+KRR, which operates in the “None” regime only.

In summary, for the chemically diverse GDB7-22-TS set, while SLATM_d allows for good performance in the “None” regime, and ChemProp in the “True” and “RXNMapper” regimes, since 3DReact can incorporate both atom-mapping information and 3D structure information, the model achieves robust performance in all three regimes, with the predicted MAEs ranging from $4.936.56\text{\,}\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{/}\mathrm{m}\mathrm{o}\mathrm{l}$.

The Cyclo-23-TS ¹¹⁷ dataset contains a single reaction class and has been previously illustrated³⁸ to show less dependence on the quality of atom-mapping than the GDB7-22-TS. For this set, 3DReact outperforms or matches the other models in all three regimes. This illustrates that a model based purely on geometry information of reactants and products, without any chemical information in the form of atom-mapping or surrogates thereof, can allow for accurate reaction property prediction. It is worth noting that atom-mapping does not improve predictions at all, i.e. there is no improvement from “None” to “RXNMapper” to “True”, even for the ChemProp model. This points to the different nature of this dataset compared to the GDB7-22-TS.

The best model is obtained with 3DReact_S in the energy mode (Figure 2d). As outlined in Section 2.2, in energy mode an energy contribution is learned for reactants’ and products’ atoms separately. In the original publication,¹¹⁷ Stuyver et al. illustrate that the activation barriers ($\Delta G^{\ddagger}$) correlate linearly with the reaction energy ($\Delta G$). Since $\Delta G$ is the difference between products’ and reactants’ energies, the energy mode is the best choice for a model learning the reaction energy, and in the case of this dataset, for $\Delta G^{\ddagger}$ too, due to its linear correlation with $\Delta G$.

Compared to SLATM_d+KRR, 3DReact in the “None” regime results in lower prediction errors for this set, illustrating that despite both models using similar information, an end-to-end model can allow for improved predictions.

The Proparg-21-TS ^{108, 69} is a small dataset for neural network standards (753 points) and therefore constitutes a challenge for the data efficiency of our model. Like the Cyclo-23-TS set, it consists of a single reaction class, i.e. enantioselective propargylation of benzaldehyde. Since the enantioselectivity is related to the barrier through an exponential relationship, it is critical to predict the barrier accurately ($\leq$ $1\text{\,}\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{/}\mathrm{m}\mathrm{o}\mathrm{l}$).⁶⁹ The “RXNMapper” regime is not available since RXNMapper cannot atom-map the reaction SMILES of this set.

In the other regimes, 3D-structure-based models lead to the best results, outperforming ChemProp by a large margin. Proparg-21-TS is particularly hard for 2D-based models³⁸ since it contains molecules of different stereochemistry but the same SMILES strings. Again trained on a single-reaction class dataset, models do not benefit from being provided the “obvious” chemical information: including true atom-maps does not decrease the error. Competing only in the “None” regime, 3DReact does not allow for a performance improvement compared to SLATM_d+KRR. Given the small size of the dataset, it is already a demonstration of data efficiency that the deep-learning model matches the prediction errors of the kernel model. Unlike for Nequip⁴³ however, the data efficiency here is not due to the equivariant molecular components (Section 3.1).

The three datasets illustrate the benefits of the flexibility of 3DReact: depending on the datasets’ particular challenges, the model exploits the available information to yield the best-performing model in almost all cases. Since the model settings (such as vector or energy mode choice) are specified as hyperparameters, the optimized version of 3DReact can emerge with minimal user intervention.

3.2.2 Extrapolative splits

Figure 5 illustrates model performance for extrapolative splits (based on scaffolds, molecular size of reactants/products, and barrier magnitude, detailed in Section 5.2). These different types of extrapolative splits are necessarily more difficult than random splits, as demonstrated by higher MAEs in Figure 5. The relative performance of the models is largely maintained in the three different extrapolation regimes compared to the interpolation regime presented in Table 2.

Figure 5: Mean absolute errors (MAEs) of predictions using three different extrapolation splits: scaffold, size-, and property-based. All datasets are compared in three atom-mapping regimes: “True”, “RXNMapper” (RXNM), and “None”, except for the Proparg-21-TS set, where RXNMapper cannot map the reaction SMILES. MAEs are averaged over 10 folds of 80/10/10 splits (training/validation/test), and error bars indicate standard deviations across folds, where applicable.

Bemis–Murcko scaffold¹²⁷ splitting clusters molecules (reactants for GDB7-22-TS and Proparg-21-TS, products for Cyclo-23-TS) based on ring systems. Test molecules may therefore appear “novel” from the point of view of the reaction graph, but will still feature distances and angles close to what the model has seen during training. Similarly for size-based splits, since there is no correlation between reactant/product size and reaction barriers, using distance information allows for stable predictions on extrapolation. Property-based splits are more challenging than the other two. For the Cyclo-23-TS and Proparg-21-TS sets, 3DReact still offers respectable errors, lower than those of the other models. For the GDB7-22-TS set however, all models result in unreasonable MAEs over $20\text{\,}\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{/}\mathrm{m}\mathrm{o}\mathrm{l}$. This points to the particular challenges of the GDB7-22-TS set and suggests an avenue for further developments of ML models for extrapolative tasks.⁸¹

Again in contrast to previous works that suggested equivariant models might be better at extrapolation tasks,^{43, 56} here we find that 3DReact offers stable extrapolation performance (particularly for size- and scaffold-based splits), but not necessarily improved extrapolation behavior compared to 2D-graph based models. This points to the different challenges in reaction property prediction. Nevertheless, Figure 5 illustrates that 3DReact is a consistently robust model for the three datasets when moving from interpolation to extrapolation regimes.

3.3 Model behavior

Since the GDB7-22-TS set has the largest chemical diversity amongst the datasets explored, studying 3DReact and baseline models SLATM_d and ChemProp on this dataset best captures the different chemical interpretation provided by these models.

Figure 6: t-SNE maps (perplexity $=64$) of the latent representations of 3DReact and ChemProp models in the “True” regime and the SLATM_d representation of the GDB7-22-TS dataset, colored by the target $\Delta E^{\ddagger}$ (upper panels) and reaction types (lower panels).

Figure 6 compares the (latent) representations of 3DReact “True”, ChemProp “True” and SLATM_d using t-SNE¹²⁸ maps. In the upper panel, we find that the quality of the correlation between the representations and the target property are aligned with the relative performance of the models (Table 2). ChemProp and 3DReact show a smooth transition of the target property, whereas the map of SLATM_d does not have a clear structure. The lower panel shows the correlation of the representations with the five most common reaction types defined by bond breaking and formation (see Section 5.4). ChemProp, as a chemically-inspired model, illustrates clear clusters in the reaction type. While SLATM_d is a geometry-based model, the binning structure used to create the representation^{8, 36} results in a clear correlation with the reaction types, since e.g. the pairwise bins naturally cluster features such as C–H bond formation or breaking. 3DReact shows the least distinct “chemical” clustering, due to the interplay of geometry and mapping information exploited in the representation.

Figure 7: Box plots illustrating how 3DReact “True” performs for the most common reaction types in the GDB7-22-TS set. 3DReact is constructed without explicit H nodes in the graphs. The boxes range from the first to the third quartile of the datapoints. The whiskers limit 90% of the datapoints and the individual points illustrate outliers. The points correspond to the test set of the first random split. The errors are given in the target standard deviation (stdev) units ($21.8\text{\,}\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{/}\mathrm{m}\mathrm{o}\mathrm{l}$).

Figure 7 shows the error distribution of predictions belonging to the same reaction classes for 3DReact “True”. 3DReact performs universally well across the different reaction types, with consistently low errors and relatively small error spread. The reactions for which the model has higher mean errors and spread ($+$H–H,$-$C–H,$-$C–H (green)) correspond to those involving C–H and H–H features. Since the model is trained without explicit H nodes in the graph, features associated with X–H bonds are included implicitly in the model. Capturing H–H bond changes will be the most challenging as these will be the least explicitly described, occurring only as initial features for neighboring nodes. Since C is the most frequently occurring element in various different configurations, capturing all the C–H features is more challenging than the O–H features for example, which will be more similar to one another. The equivalent plot for the model trained with explicit H nodes is shown in Figure LABEL:S-fig:box_plot_with_H, illustrating that the error spread reduces for the reaction types involving C–H and H–H features. Note that 3DReact without explicit Hs still leads to performance comparable to the variant with explicit Hs (Section LABEL:S-sec:hydrogens).

3.4 Geometry quality

In order to illustrate that 3DReact does not require high-quality molecular structures to be used in an out-of-sample scenario, we train and test a model using lower-quality GFN2-xTB¹²⁹ (xTB) geometries to predict higher-level barriers (CCSD(T)-F12a/cc-pVDZ-F12//$\omega$B97X-D3/def2-TZVP for GDB7-22-TS, B3LYP-D3(BJ)/def2-TZVP//B3LYP-D3(BJ)/def2-SVP for Cyclo-23-TS and B97D/TZV(2p,2d) for Proparg-21-TS). The results are illustrated in Figure 8 for the three datasets with DFT and xTB geometries, and compared to the SLATM_d+KRR model in the same settings. 3DReact benefits from a lower sensitivity to the geometry quality compared to the pre-designed representation SLATM_d combined with KRR, across the three datasets.

Figure 8: Mean absolute errors (MAEs) for predictions using either the provided geometries ($\omega$B97X-D3/def2-TZVP for GDB7-22-TS, B3LYP-D3(BJ)/def2-SVP for Cyclo-23-TS, B97D/TZV(2p,2d) for Proparg-21-TS) (DFT) or lower-quality GFN2-xTB (xTB) geometries. MAEs are averaged over 10 folds of random 80/10/10 splits (training/validation/test), error bars showing standard deviations across folds. Note that for GDB7-22-TS and Cyclo-23-TS datasets the DFT results are different from those presented in Section 3.2.1 because here they are obtained on the same subset as the xTB results (see Section 5.4).

For the GDB7-22-TS set, there is a negligible difference in model performance moving from DFT to xTB geometries. The xTB geometries are a good proxy for the DFT ones here, since this set consists of small, charge-neutral organic molecules, which are largely well-described by semi-empirical methods. For the Cyclo-23-TS set, while the molecules are still organic, they are larger than those in the GDB7-22-TS set, and there is a greater divergence between the GFN2-xTB and DFT geometries, resulting in a larger deterioration with these structures. Figure LABEL:S-fig:cyclo_rmsd demonstrates that when using the model trained on xTB geometries, barrier predictions for molecules with poorer geometries (i.e., higher RMSD of xTB vs. DFT geometries) are not necessarily worse than those on molecules with better geometries. Instead, there is a consistent decline in model performance when training with xTB geometries and predicting DFT barriers.

The Proparg-21-TS set is the most complex of the three for GFN2-xTB, since these systems with charged organosilicon compounds differ considerably from those used to parameterize semi-empirical methods or force fields. As described in Section 5.4, unlike for the other datasets where we generate an initial structure from SMILES using force fields, for this set it is impossible and we instead generate xTB geometries from the DFT ones. While this is not a feasible geometry generation pipeline for out-of-sample predictions, it still demonstrates how different methods perform with high and low-quality geometries. Here, we see that 3DReact is less sensitive than SLATM_d+KRR and the variation trained with lower quality geometries still offers competitive errors ($0.48\pm 0.05$ kcal/mol for the “None” model).

4 Conclusions

The accurate and reliable prediction of reaction barriers across diverse sets of chemical reactions remains an open challenge in computational chemistry. We contribute to this domain by introducing 3DReact, a geometric deep learning model constructed from the 3D coordinates of reactants and products. We show that the invariant model (vs. the equivariant version) is already sufficient for currently available reaction datasets. Existing models ChemProp and SLATM_d+KRR exhibit impressive performance for atom-mapped, chemically diverse datasets and stereochemistry-sensitive datasets, respectively. 3DReact offers a hybrid model that can optionally incorporate mapping information alongside geometries, enabling robust performance across different dataset types and atom-mapping regimes. 3DReact also allows for a reduced sensitivity to the training geometry quality (i.e., xTB vs. DFT level) compared to SLATM_d+KRR. Predictions are stable both when moving to molecular size- or scaffold-based splits. Altogether, 3DReact presents a flexible framework for accurate prediction of activation barriers across chemical reaction datasets. Despite the proposed developments, challenges remain for ML predictions of energy barriers, particularly in integrating them within experimental settings. This work is a step toward their reliable application.

5 Methods

5.1 Datasets

We test 3DReact on three datasets of reaction barriers previously used to benchmark reaction representations.³⁸ The term “reaction barrier”, used interchangeably with “activation energy” and “activation barrier” is the energy difference between the energy of the optimized TS and the optimized reactants. Note that depending on the dataset, some provide purely electronic energies (labelled $\Delta E^{\ddagger}$) and others — Gibbs free energies (labelled $\Delta G^{\ddagger}$). In all datasets, optimized three-dimensional structures of reactants and products are provided, which are used to train models and make predictions. The activation barrier is not a direct function of these structures, but using the TS structure to make predictions removes the utility of the ML models vs. direct computation of the TS. Thus we use an implicit interpolation of reactants’ and products’ structures as a proxy for the TS as in previous works.^{36, 69, 38}

The GDB7-22-TS ¹¹⁶ dataset consists of close to $12\,000$ diverse organic reactions automatically constructed from the GDB7 dataset^{130, 131, 132} using the growing string method¹³³ along with corresponding energy barriers ($\Delta E^{\ddagger}$) computed at the CCSD(T)-F12a/cc-pVDZ-F12//$\omega$B97X-D3/def2-TZVP level. The dataset provides atom-mapped SMILES, with “True” maps derived from the transition state. For $43$ reactions out of $11\,926$, one of the products’ SMILES represents a molecule different from the xyz structure. These reactions were therefore excluded from the dataset, leading to a modified GDB7-22-TS set used here.

While there are no pre-defined classes for all the reactions in the GDB7-20-TS¹²³ or GDB7-22-TS ¹¹⁶ sets, Grambow et al. ⁷⁰ split the dataset into reactions undergoing certain bond changes: for example, the most common type was breaking of a C–H bond ($-$C–H) and a C–C bond ($-$C–C) in the reactants and formation of a C–H bond ($+$C–H) in the products, giving the reaction type signature $+$C–H,$-$C–C,$-$C–H. Here, we extract similar reaction types by comparing the connectivity matrices from atom-mapped reaction SMILES of reactants and products (ignoring bond orders). The most abundant reaction types in the dataset are $+$C–H,$-$C–C,$-$C–H (1667 reactions), $+$H–N,$-$C–H (633), $+$C–H,$-$C–H (619), $+$H–O,$-$C–H,$-$C–O (599) and $+$H–H,$-$C–H,$-$C–H (517).

The original Cyclo-23-TS ¹¹⁷ dataset encompasses $5\,269$ profiles for $[3+2]$ cycloaddition reactions with activation free energies ($\Delta G^{\ddagger}$) computed at the B3LYP-D3(BJ)/def2-TZVP//B3LYP-D3(BJ)/def2-SVP level in water using the SMD continuum solvation model. The dataset provides atom-mapped SMILES with “True” maps for heavy atoms derived from either the transition state structure or heuristic rules. For the regime with explicit hydrogen atoms, we atom-mapped the xyz files by matching the reactants, given in two separate files, to the provided transition state structure, which closely resembles the two reactants and has the same atom order as in the product. This was done with a labelled graph matching algorithm as implemented in NetworkX.^{134, 135} The algorithm is unaware of chirality, double-bond stereochemistry or conformations, and thus may lead to not exactly correct atom-mappings. We also found that in four reactions, the product SMILES and xyz files depict different species, thus the set was reduced to $5\,265$ reactions.

The Proparg-21-TS dataset^{108, 69} contains 753 structures of intermediates before and after the enantioselective transition state of benzaldehyde propargylation, with activation energies ($\Delta E^{\ddagger}$) computed at the B97D/TZV(2p,2d) level. SMILES strings (“fragment-based” SMILES) and “True” atom-maps are not provided with the original dataset, these are taken from Ref. 38.

RXNMapper¹⁰⁷-mapped versions of GDB7-22-TS and Cyclo-23-TS were obtained with the python package rxnmapper (version 0.3.0), using the default settings. The Proparg-21-TS set cannot be mapped, because the underlying libraries cannot process its SMILES string.³⁸ Since RXNMapper sorts molecules in case of multiple reactants and/or products, which would complicate SMILES–xyz matching (see Section 5.3 below), we used a locally modified version that does not change the molecule order (the patch file is provided in the project repository at https://github.com/lcmd-epfl/EquiReact/tree/9d78892fe/data-curation/rxnmapper).

5.2 Data splits

For each dataset and splitting type, identical data splits were used for all the models compared. In each case, ten different splits are constructed with different random seeds.

Three different types of extrapolation split were used: scaffold-, molecular size- and property-based. Scaffold splitting^{136, 137} clusters molecules based on their 2D backbones (such as Bemis–Murcko scaffolds¹²⁷) and ensures that the clusters (scaffolds) belonging to the training, validation, and test sets do not overlap. Size-based splitting organizes the splits such that the reactions of the smallest molecules are in the training set and the reactions of the largest molecules are in validation and test. With property-based splits, one trains on reactions with higher barriers and predicts on reactions with lower barriers. This choice of splits reflects the relevant out-of-sample cases: larger molecules are more expensive to compute, and reactions with smaller barriers are desirable. Size- and property-based splits can also be organized in reverse order, where larger molecules are in the train set and smaller in test, or reactions with lower barrier in train and higher barrier in test.

For molecular size- and scaffold-based splits, the initial data shuffling affects the composition of the datasets. The non-zero standard deviations for property-based splits with neural networks arise from different organization of the datapoints into batches.

5.3 Matching SMILES strings to xyz geometries

3DReact makes use of both the graph structure of a molecule (as provided in the SMILES string) and the three-dimensional structure (in the xyz). The atoms in the graph are associated with the atomic coordinates provided in the xyz file. Thanks to the way the GDB7-22-TS dataset¹¹⁶ was generated, the atomic coordinates can be easily matched to SMILES which in turn allows to atom-map reactants to products. However, we also tested RXNMapper-mapped SMILES which do not respect the same constraints. Therefore, for consistency, we use a SMILES–xyz matching procedure detailed below.

We construct molecular graphs from xyz using covalent radii and matched them to RDKit¹²⁰ molecular graphs obtained from SMILES with a labelled graph matching algorithm as implemented in NetworkX.^{134, 135} This procedure is however unaware of chirality and double-bond stereochemistry, thus some of the matches might be incorrect. Still, it provides a flexible method that can be applied to any dataset consisting of SMILES strings and xyz files.

The same procedure was applied to the Cyclo-23-TS dataset in the few cases when the canonical SMILES have a different atom ordering than xyz.

5.4 xTB geometry generation

For the GDB7-22-TS and Cyclo-23-TS datasets, the starting structures were generated from SMILES using the distance-geometry embedding implemented in RDKit¹²⁰ with the srETKDGv3 settings.¹³⁸ Ten conformations were produced per molecule, which were then energy-ranked with the MMFF94 implementation¹³⁹ in RDKit, defaulting to UFF in case of missing parameters. The lowest energy conformer was retained. For the Proparg-21-TS set, the original B97D/TZV(2p,2d) geometries were used as a starting point, because the stereochemical and conformational diversity of this set cannot be completely encoded with SMILES. Therefore MMFF94 will fail to generate an initial geometry from SMILES.

For all the sets, the starting structures were optimized at the GFN2-xTB semi-empirical level of theory¹²⁹ at the “loose” convergence level for a maximum of 1000 iterations using xTB¹⁴⁰ version 6.2 RC2. For $969$ reactions of the GDB7-22-TS set and $491$ reactions of the Cyclo-23-TS set, at least one of the participating molecules either could not converge to any reasonable configuration or converged to a structure not matching the SMILES. These reactions were excluded from the geometry quality tests (Section 3.4).

5.5 Model training

3DReact was trained using the Adam optimizer ¹⁴¹ with initial learning rate and weight decay parameters as hyperparameters. The learning rate was reduced by 40% after $60$ epochs of no improvement in the validation MAE, as in Ref. 121. Models were trained for max. $512$ epochs, using early stopping after $150$ epochs of no improvement. The model with the best validation score was then used to make predictions on the test set.

The optimal model hyperparameters were searched within the following values: learning rate $\in[$5\text{$\cdot$}{10}^{-5}$,${10}^{-4}$,$5\text{$\cdot$}{10}^{-4}$,${10}^{-3}$]$; weight decay parameter $\in[${10}^{-5}$,${10}^{-4}$,${10}^{-3}$,0]$; node and edge features embedding size $n_{s}\in[16,32,48,64]$; $\ell{=}1$ hidden space size $n_{v}\in[16,32,48,64]$; number of edge features $n_{g}\in[16,32,48,64]$; number of convolutional layers $n_{\mathrm{conv}}\in[2,3]$; radial cutoff $r_{\max}\in[2.5,5.0,10.0]$; maximum number of atom neighbors $n_{\mathrm{neigh}}\in[10,25,50]$; dropout probability $p_{d}\in[0.0,0.05,0.1]$; sum_mode $\in$ [node, both]; combine_mode $\in$ [mlp, diff, mean, sum]; graph_mode $\in$ [energy, vector].

The hyperparameter search was done for the equivariant model EquiReact_S (without attention or mapping) using Bayesian search as implemented in Weights & Biases.¹⁴² Hydrogen atoms were excluded from the graphs. Sweeps were run for 128 epochs for the GDB7-22-TS and Proparg-21-TS sets, and for 256 epochs for the Cyclo-23-TS set on the first random split. The parameters resulting in the best validation error, summarized in Table LABEL:S-tab:model-params, were used for all the other model settings.

5.6 Baseline models

The ChemProp model¹⁰³ is based on a CGR built from atom-mapped SMILES strings of reactants and products, which is then passed through the directed message-passing neural network chemprop^{137, 71, 103} (version 1.5.0). The hyperparameters are taken from Ref. 38.

Molecular SLATM vectors were generated using the qml python package¹⁴³ before being combined to form the reaction version SLATM_d. SLATM_d is used with kernel ridge regression (KRR) models. The kernel functions and widths, and regularization parameters, were optimized on the first of the ten random splits, in line with how the hyperparameters were optimized for 3DReact. Unlike 3DReact, the hyperparameters for DFT and xTB geometries were optimized separately.

Data and Software Availability statement

The code is available as a GitHub repository at https://github.com/lcmd-epfl/EquiReact. The versions of the datasets used, as well as any processing applied to them, can be found in the same repository. The unprocessed results are available in the same same repository as well as at https://wandb.ai/equireact.

Supplementary Information is provided in the freely available file equireact_si.pdf, detailing the architecture of the molecular channels (Section LABEL:S-sec:molecular_channels), the 3DReact hyperparameters (Section LABEL:S-sec:model-params), the RMSE analogue of Table 2 (Section LABEL:S-sec:rmse), the discussion of the model with a cross-attention surrogate for atom-mapping (Section LABEL:S-sec:cross), extrapolation studies (Section LABEL:S-sec:extrapolation), some illustrative correlation plots for the GDB7-22-TS set (Section LABEL:S-sec:gdb_outliers_corr), the model performance with and without explicit hydrogen atoms (Section LABEL:S-sec:hydrogens), and the geometry sensitivity analysis for the Cyclo-23-TS set (Section LABEL:S-sec:geom_cyclo).

Author Information

Author contributions

P.v.G., K.R.B., and C.B. conceptualized the project. 3DReact and support codes were written and run by K.R.B. and P.v.G., with design suggestions from C.B., V.R.S., and R.L. Results were analyzed by P.v.G., K.R.B., V.R.S., R.L., and C.B. xTB computations were run by R.L. The original draft was written by P.v.G. and K.R.B. with reviews and edits from all authors. C.C. and A.K. provided supervision and acquired funding.

Conflict of interest

The authors have no conflicts to disclose.

The authors thank Liam Marsh and Yannick Calvino Alonso for helpful discussion and comments on the text. P.v.G., C.B., V.R.S., R.L., A.K., and C.C. acknowledge the National Centre of Competence in Research (NCCR) “Sustainable chemical process through catalysis (Catalysis)”, grant number 180544, of the Swiss National Science Foundation (SNSF) for financial support. K.R.B. and C.C. were supported by the European Research Council (grant number 817977) and by the National Centre of Competence in Research (NCCR) “Materials’ Revolution: Computational Design and Discovery of Novel Materials (MARVEL)”, grant number 205602, of the Swiss National Science Foundation.