Thousands of reactants and transition states for competing E2 and S$_\mathrm{N}$2 reactions

In this study we used delta machine learning (Δ-ML) in kernel ridge regression (KRR) implemented in the $\mathrm{QMLcode}$ [58]. Kernel based methods were introduced in the 1950 s by Kriging et al [59]. KRR uses as input a kernel function with the feature vector $\mathbf{x}$ to learn a mapping function to a property $y_{q}^{\mathrm{est}}(\mathbf{x}_{\,j})$ given a training set of N reference pairs $\{\mathbf{x}_{\,i}, y_i\}^{N}$:

$$\begin{equation} y_{q}^{\,\mathrm{est}}(\mathbf{x}_{\,j}) = \sum^{N}_{i} \alpha_i k(\mathbf{x}_{\,i}, \mathbf{x}_{\,j}) \end{equation} \tag{ 1 }$$

where α is the regularization coefficient and $k(\mathbf{x}_{\,i}, \mathbf{x}_{\,j})$ a Gaussian kernel element:

$$\begin{equation} k(\mathbf{x}_{\,i}, \mathbf{x}_{\,j}) = \mathrm{{exp}}\left( -\frac{||\mathbf{x}_{\,i} - \mathbf{x}_{\,j}||_2^{2}}{2\sigma^{2}}\right) \end{equation} \tag{ 2 }$$

A more detailed discussion of the KRR method employed in this work and pertinent references can be found in Heinen et al [60]. In the context of Δ-ML, the procedure stays the same and only the property (y) changes from a molecular property to a difference in properties, e.g. from $y^{\mathrm{est}} \mathrel{\widehat{ = }} E_{\mathrm{a}}$ to $y^{\mathrm{est}} \mathrel{\widehat{ = }} \Delta E_{\mathrm{a}}$.

The feature vector or representation $\mathbf{x}$ we used is one-hot encoding [61], which is a bit vector. For every substitution site Rk, nucleophile Y and leaving group X, we denote presence of a given combination with ones. In our case, this means that for any transition state, six out of the 27 entries of the representation vector are ones, the rest zeros.

We started from the unsubstituted case fluoroethane optimized with openbabel [62] using the universal force field (UFF) [63] and functionalized the substituent sites Rk in figure 1 using the C++ interface of openbabel. Again, each resulting structure was optimized with UFF to remove potential bad contacts. Using the Experimental-Torsion Knowledge Distance Geometry (ETKDG) method as implemented in RDKit [64], we searched for 1,000 conformer geometries. They subsequently were ordered by UFF energy. Starting from the most stable conformer, all those configurations were included in the followings steps if and only if their root mean squared difference (RMSD) to the previously accepted configuration was at least 0.01 Å or the energy difference between the two was at least 0.1 kcal mol⁻¹.

The resulting conformer candidate configurations were relaxed at MP2/6-311G(d) level with ORCA 4.0.1 [48–50, 65–67] to be compatible with the level of theory to be employed for the transition state search. For each of these minimized configurations, all possible nucleophiles given in table 1 were placed along the expected axis of the CH bond in figure 3. With the nucleophile being constrained to that axis, the geometries were optimized to obtain an estimate of the reactant complex geometry.

For each of these reactant complexes, we subsequently lifted the constraint and relaxed further. This was helpful as the potential energy landscape around the reactant complex is comparably shallow and therefore direct optimization to the free reactant complex was often ineffective.

Each unconstrained reactant complex was validated using a variety of geometrical criteria to ensure that the more than 100,000 minimum energy geometries represented meaningful configurations. The overall procedure is shown in figure 2. First, we required the reactant complex to constitute two fragments based on the topology obtained from MDAnalysis [68] where one fragment needed to be of exactly one atom, i.e. the nucleophile. This is to avoid erroneous fragmentation where e.g. a proton is abstracted from the reactant. In the case of E2 reactions, we required that the angle C–H$\cdots$Y must not be smaller than 178 degrees since configurations with larger angles indicate trapping of the nucleophile by other hydrogen atoms of the reactant not involved in this particular reaction channel.

Zoom In Zoom Out Reset image size

For S$_\mathrm{N}$2, more validations are required.

• The C$\cdots$Y distance had to be at least 1.14 Å, 1.41 Å, 1.86 Å, and 2.04 Å for hydrogen, fluorine, chlorine, and bromine, respectively. This avoids configurations that are actually product complexes. Due to the low activation energy for many such cases, a geometry optimization can end up in a product complex minimum from a reactant complex initial guess.
• To avoid trapping of the nucleophile by reactant hydrogen atoms, the distance between the nucleophile and the closest hydrogen of the reactant is required to be at least 0.78 Å, 0.96 Å, 1.33 Å, and 2.48 Å for hydrogen, fluorine, chlorine, and bromine, respectively.
• Since the S$_\mathrm{N}$2 reaction requires nearly planar bonds for the reaction center, we require that the angle X$\cdots$C$\cdots$Y must be at least 178 degrees.
• We avoid artificially stretched geometries by requiring no carbon-carbon distance to be within 1.65–2 Å and no nitrogen-oxygen distance to be within 1.5–2.5 Å.

Whenever these validation steps were successful, the lowest such minimum from all conformers investigated is considered to represent the reactant complex. Otherwise, the lowest energy configuration from the constrained optimization is taken as an approximation of the reactant complex. In the latter case, Δ-ML [36] was employed to estimate the residual relaxation energy between the constrained and unconstrained reactant complex.

Duplicate reactant and product conformer geometries were identified using the FCHL19 [26] representation. By that measure, only unique geometries are retained. This test was not applied to reactant complexes as their local minima energies and geometries can be very similar yet distinct.

Using Gaussian09 [69] for an initial transition state geometry with B3LYP/6-31G* [70–75] and subsequently ORCA 4.0.1 for a final transition state with MP2, we first found the transition state of the unsubstituted case with chloride as nucleophile. Functionalization followed the same procedure as for the reactants. Using these starting geometries, transition states were obtained via eigen mode following as implemented in ORCA. After a transition state was found, the local Hessian matrix was obtained from a numerical frequency calculation by finite displacements as implemented in ORCA.

Once a transition state was found for a combination of the four substituents, this geometry was employed as starting geometry for further transition state searches for missing cases where exactly one out of the four substituents was different from the case where a validated transition state has been found. This scheme was used only for those molecules where the substituent that was to be replaced did not have the same functional group as the neighbouring substituent on the same carbon atom. For some cases, this procedure was employed several times in a row, each time resulting in an additional set of transition states which served as starting guesses. Similarly, the nucleophile of validated transition states was replaced to obtain promising starting geometries for the transition state search.

Once the transition state geometry has been found for any potential reaction target, the Hessian was evaluated to ascertain that the geometry in fact is a transition state with exactly one imaginary frequency. We only included a transition state in our data set if this frequency was at least 400 cm⁻¹ and that the resulting motion corresponding with this one normal mode was as shown in figure 3 (left column). The ethane skeleton features two carbon atoms Ck, where the one with substituents R1 and R2 is numbered C1. For the E2 transition state, X, Y and the hydrogen atom were displaced along the normal mode and checked if the distances C2-H as well as C1-X were larger and the C2-Y distance was smaller compared to the non-displaced geometry. In the S$_\mathrm{N}$2 transition state, the nucleophile and leaving group were displaced along the normal mode and C1-X was compared to C1-Y.

Zoom In Zoom Out Reset image size

While the investigation of the normal modes alone ensures that the vibrational motion belongs to the main configurational change the molecule undergoes during each reaction, it is not a sufficient criterion that this particular transition state geometry actually connects reactant and product. We use the intrinsic reaction coordinate (IRC) [76] as final criterion to ensure that the transition state indeed connects a valid reactant complex with a valid product for the reaction in question. The IRC is commonly employed to find a reaction pathway starting from a transition state. The Cartesian IRC is given by the steepest descent path in forward and backward direction of the reaction. We use steepest descent as implemented in ORCA to trace the Cartesian IRC. If the energy curvature near the transition state and along the reaction coordinate is small, steepest descent paths can become subject to numerical instabilities. To avoid this issue, we approximate the IRC close to the transition state by a line scan in either direction based on the normal mode displacement of imaginary frequency. From the final point of the line scan, a regular steepest descent is followed until a local minimum has been reached.

Since the sign of the normal mode of imaginary frequency is not fixed with respect to the direction of the reaction, we analyze the minimum energy endpoints of the IRC to classify them as either close to reactant or close to product based on the bond length as shown in figure 3. If and only if exactly one of the endpoints is found to be close in geometry to a reactant configuration and the other is found to be close in geometry to the product configuration, the corresponding transition state is included in our data set. To test whether the configurations are close to reactant or product, we measured C2-H distances for the E2 case and C1-X and C1-Y distances for the S$_\mathrm{N}$2 reaction to ensure bonds have been broken as shown in figure 3.

For cases where several validated transition states for the same reaction have been found, we consider the lowest one for the reaction barrier.

Finally, we performed single-point DF-LCCSD/cc-pVTZ calculations, as implemented in Molpro2018 [77–83] using the extremal geometries as obtained with MP2/6-311G(d). All in all, the complete generation of the data set took about 2.8 million core hours.

Our resulting data set contains 4,466 validated transition state geometries, of which 2,785 are for S$_\mathrm{N}$2 (TS$_\mathrm{S}$) and 1,681 for E2 (TS$_\mathrm{E}$). Based on 26,997 reactant conformers (R$_\mathrm{S,E}$), we identified 81,950 constrained reactant complexes for E2 (R'$_\mathrm{E}$) and 57,642 constrained reactant complexes for S$_\mathrm{N}$2 (R'$_\mathrm{S}$) which in turn have been refined to yield 2,030 unconstrained reactant complexes for E2 (R'$_\mathrm{E}$) and 1,532 unconstrained reactant complexes for S$_\mathrm{N}$2 (R'$_\mathrm{S}$). Finally, we have found 15,706 S$_\mathrm{N}$2 product conformers (P$_\mathrm{S}$) and 9,588 E2 product conformers (P$_\mathrm{E}$). All geometries are calculated at MP2/6-311G(d) level of theory and given as XYZ files in this work. Two additional files specify all individual energies and activation energies, respectively. The labels in the text files relate to the labels in table 1.

All data is available in the materials cloud (https://doi.org/10.24435/materialscloud:sf-tz).

As shown in figure 4, we were able to find many transition states for a variety of substituents, nucleophiles and leaving groups. This means that we have reached a substantial coverage of the chemical space in question, which is key for machine learning. The challenge here is the low success rate of the transition state search which might have been the key reason why such data sets have not yet been published earlier. In particular, machine learning models will benefit from the comparably low noise in the data set coming from our validation procedure. Moreover, the data set features many different combinations of substituents such that there is considerable promise that their interplay for the competing reactions can be analysed and understood.

Zoom In Zoom Out Reset image size

As in any iterative optimization scheme, convergence thresholds influence the final results. This is the case for a transition state search as well and might potentially give rise to some small noise in the transition state geometries. Since we calculated the explicit Hessian matrix, we know that the transition state geometries reported in this data set are indeed saddle points, and that their mass-weighted normal modes represent the concerted rearrangement expected for E2 and S$_\mathrm{N}$2 reactions. Together with the tight convergence criteria required for transition state optimization, this means that our data set contains only highly compatible transition state geometries for all the validated combinations of substituents, nucleophiles and leaving groups.

This is demonstrated in figure 5 which shows a scatter plot of the most important internal coordinates for the transition states. The reduction of dimensionality from the more complex 3D geometry is obtained by placing one of the two central carbon atoms in the origin and then aligning the carbon-carbon bond along one Cartesian axis. The other markers then show the position of one atom for each transition state found. For E2 reactions, the transition state geometry has been rotated such that all three points shown in the corresponding panels are exactly within one plane. For the S$_\mathrm{N}$2 case, this is not possible, as the four atoms in question are not necessarily exactly in a plane even though they are very close to that. For this panel, the projection on the fitted plane through all four points is shown. For the internal carbon-carbon bond, the variance of the bond length is significantly higher for E2 than for S$_\mathrm{N}$2, as shown in figure 5. This can be explained by the nature of the two reactions: While E2 consists of a concerted action on both carbons, S$_\mathrm{N}$2 happens only at one of the two carbons. We also see that each element for the nucleophile and leaving group has its own distribution of positions relative to the two central carbon atoms. This distribution reflects the impact of the different substituents on each transition state geometry. It is interesting that fluorine atoms exhibit much less spatial variation as leaving group than other halogens for E2 while this is not at all the case for the role of fluorine as nucleophile in the very same reaction. This is likely attributed to the comparably short bond distance of fluorine for the leaving group, since in the case of the nucleophile this distance is increased due to one intermediate hydrogen atom between the central carbon and the nucleophile. The reduced distance in the former case then would lead to a more pronounced Coulombic interaction with the molecule, effectively restraining the fluorine atom to a smaller volume of configurational space.

Zoom In Zoom Out Reset image size

The centers of the positional distributions of the three halogens as leaving group increase with the period of the element, which is in line with typical bond radii for these elements. This is more pronounced in the case of the nucleophiles in E2 reactions where the intermediate hydrogen atom reduces the interaction between nucleophile and molecule. The result is that the nucleophile positions are spread out on arcs around the central carbon with most of the positional freedom captured by the intermediate hydrogen atom. Again, the radii of the halogen arcs follow the period of the elements, while a hydrogen as nucleophile is most flexible in regards to its distance from the central carbon.

For the distribution of internal coordinates for the S$_\mathrm{N}$2 reaction in figure 5, two features are most striking: the triangular domain of the positions of halogenic nucleophiles and the bimodal distribution of hydrogens in the same case which in turn is mirrored in a bimodal distribution for the leaving group positions for all elements.

The triangular domain for halogenic nucleophiles in figure 5 can be explained by their electrostatic interaction with the reactant molecule in the gas phase. For the transition state to be a saddle point, all but one degrees of freedom must yield an increase of energy. At the tip of the triangular domain, there are three bounds to observe. First, if the distance to the carbon forming the reaction center would decrease, then the binding energy gain would become dominant, so this distance needs to be slightly above the equilibrium bond length. Secondly, the direction towards the planar substituents R1 and R2 would reduce the distance between the partially negatively charged nucleophiles and the partially positively charged hydrogen atoms of the substituents. This Coulombic interaction is more pronounced in gas-phase and restricts the possible geometries for transition states in this direction. Finally, pushing more towards the other carbon atom of the reactant skeleton (upwards in figure 5), would be unfavourable in the sp² hybridisation of the reaction center. Only for larger distances of the nucleophiles to the reaction center, deviations from the last two constraining factors become possible, hence the triangular shape of the domain for each element.

The bimodal distribution of the hydrogen nucleophiles for S$_\mathrm{N}$2 as shown in figure 5 correlates with the leaving group in the corresponding reaction. Only leaving groups of chlorine and bromine allow a distance C2-H larger than 2 Å. This could be linked to the substantially higher electronegativity of fluorine, pulling more of the C2 electron cloud towards the leaving group, allowing for a shorter distance to the H⁻.

Results such as the triangular domains and the bimodal distribution can be easily identified in large homogeneous data sets such as this one and can be interesting test cases for machine learning models for phenomena resulting from the complex interplay of competing physical interactions.

Based on the conformational search for the reactant geometries and the validated transition states, we could calculate activation energies for both reactions. Figure 6 shows the broad distribution of said activation energies which span about 50 kcal/mol. In general, E2 activation energies are lower than S$_\mathrm{N}$2 activation energies. Since the activation energies are defined as the difference in energy between the transition state and the reactant complex, the nature of the reactant complex is highly relevant. This is exemplified by the significant portion of negative activation energies if we consider the constrained approximation of the reactant complex alone (panel (a) in figure 6).

Zoom In Zoom Out Reset image size

These spurious negative activation energies result from two aspects: the finite number of conformers tested as potential reactant complex geometry and the constraint enforcing the characteristic alignment of the nucleophiles with the molecule when forming a reactant complex. To alleviate the impact of the former effect, we searched for more conformers until the number of negative activation energies could not be reduced any further despite testing of additional conformers. Here, the small size of the molecular skeleton was helpful, as only a few conformers can be realized for each molecule in our chemical space. We dealt with the second reason for negative activation energies by removing the constraint for the characteristic alignment of the nucleophiles with the molecule. This constraint was needed initially to ensure that the relaxation (described in the Methods section above) did not converge to an irrelevant reaction complex where the nucleophiles would be trapped by the partially positively charged hydrogen atoms of the substituents. Since the minimum of the reaction complex is only shallow, this initial constraint drastically improved the success rate of finding reactant complexes matching the reaction mechanism.

Relaxing the reactant complexes further without the constraint again bears the risk of the substituents trapping the nucleophiles. Consequently, many but not all reactant complexes could be refined this way: 301 and 348 targets for E2 and S$_\mathrm{N}$2, respectively. We expect that turning the constraint into a restraint that subsequently is reduced during the minimization until the unconstrained minimum is found could be one route to identify the correct relaxation energy for all reactant complexes in our data set. However, this would be extremely costly and is subject to many degrees of freedom, like the speed at which the restraint is removed such that this route is not feasible for the thousands of reactant complexes we have in our data set. Therefore, we trained a one-hot-encoding KRR machine learning model to take the explicit relaxation energies we have found and to predict the relaxation energies for the remaining compounds. These relaxation energies span about 15 kcal/mol. We could machine learn the relaxation energy down to prediction errors of 1.5 and 1.8 kcal/mol (for 280 randomly chosen training instances) for two separate models for S$_\mathrm{N}$2 and E2 reactions, respectively (see inset panel (b) of figure 6). This is much less than the expected error of the quantum chemistry method that we use, MP2. We do expect that more sophisticated machine learning methods could possibly improve upon this accuracy.

Panel (b) in figure 6 shows the activation energies for those barriers where we were able to find the explicit minimum geometry for the unconstrained reactant complex. The fact that this exhibits nearly no negative activation energy is in line with our observation that searching for additional conformers as basis for the reactant complex did not yield any further change to the activation energies. Using the explicitly calculated relaxed reactant complexes where available and including a machine learned relaxation energy in the activation energy for all other reactions, we obtain our final MP2/6-311G(d) and ML corrected MP2/6-311G(d) numbers for the activation energy, shown in panel (c) of figure 6 which now span 60 kcal/mol for S$_\mathrm{N}$2 reactions and 50 kcal/mol for E2 reactions.

Comparing panels a) and c) in figure 6 shows how the number of negative activation energies has been greatly reduced by removing the constraint on the reactant complex, confirming that this was the main reason for negative activation energies in the initial case of panel (a). Calculating activation energies directly from the reactants at infinite separation is no substitute for this complicated procedure of correcting for the constraint impact. This is due to the significant interaction energy gained in forming the reactant complex in the E2 case where the negatively charged nucleophile approaches a hydrogen. For S$_\mathrm{N}$2, in few cases the reactant complex might be higher in energy than the reactants at infinite separation.

Given the documented quality of MP2 geometries for substitution reactions [8], the main difference to higher level of theory than MP2 is expected to come from higher-quality energies for MP2 geometries. Since higher level of theory calculations are not affordable in the context of the geometry optimizations for this many configurations, additional single points on top of MP2 geometries recover at least a substantial part of the difference in the potential energy landscape. For those cases where we have both the transition state and the unconstrained reactant complex, we performed DF-LCCSD/cc-pVTZ calculations. The explicit data is shown in panel (d) of figure 6. The difference to the MP2 data however is more interesting and shown in panel (e) of the same figure. While the distribution of the corrections is centered around zero, the typical correction is on the order of a few kcal/mol with only few substantially larger values.

Explicit calculations of the LCCSD energy are only accessible for cases where we have an explicit molecular geometry. If the unconstrained geometry optimization did not successfully find the shallow minimum of the reactant complex, then this explicit molecular geometry is not available. To extend the coverage of the LCCSD correction which improves the accuracy of the activation energy data, we built a one-hot encoding machine learning model that predicts the LCCSD energy for the missing geometries. This Δ-ML approach exhibits learning with an error of less than 1 kcal/mol (S$_\mathrm{N}$2) and 1.5 kcal/mol (E2) after training on 280 instances. After this second step, we obtain our final activation energies which have a slightly broader distribution than before, shown in panel (f) of figure 6. It is interesting to note that this final activation energy distribution of the E2 is dramatically more skewed towards very small values than the S$_\mathrm{N}$2 which appears to be more normally distributed. This could be due to the symmetry in the case of the S$_\mathrm{N}$2 (as also shown in figure 5) where one covalent bond is broken as the other is formed. The E2 reaction is less symmetric, effectively breaking one single bond while forming a double bond. The structural lack of symmetry is also on display in figure 5.

We also note that the learning curves for the activation energy of E2 display a higher off-set than for S$_\mathrm{N}$2 even though, the E2 data has a smaller magnitude and variance. This latter aspect could be due to some extreme outliers in the E2 data set for which values larger than 50 kcal/mol have been observed, introducing severe bias in the mean absolute error. A median error measure might be better tempered for such a data set.

Panel (f) of figure 6 shows some remaining negative activation energies. For S$_\mathrm{N}$2, there are 43 such negative energies, all but one of which are from machine learning predictions only. For E2, there are 120 such negative energies in total, 79 of which come from machine learning predictions. Therefore, for the majority of cases the machine learning model needs improvement, possibly by adding more explicit unconstrained reactant complexes. The cases where the explicitly calculated activation energies are still negative likely come from a finite search of conformer geometries, meaning that some unconstrained reactant complex minima have not been found. In our data set, we include these negative activation energies such that future machine learning models correcting e.g. the constrained to unconstrained relaxation can test whether they improve upon our approach.