Relative rates of addition of carbanions to substituted nitroarenes: can quantum chemical calculations give meaningful predictions?

Abstract

Computational description and kinetic properties based on density functional theory methods of the key step of the addition reaction between a model nucleophile and nitroaromatic ring in positions occupied by hydrogen are presented. A wide series of DFT functionals (PBE0, B3LYP, ωB97XD, M062X, PBE1PBE-D3, B3LYP-D3 and APFD) was used to track the influence of functional groups in nitroaromatic rings on reaction activation barriers. The comparison of experimentally determined relative rates of nucleophilic addition and their calculated thermodynamic counterparts at various positions in a series of ortho-, meta- and para-substituted nitroarenes are provided. It was shown that different DFT methods provide a good correlation between computed thermodynamic parameters and logarithms of experimental relative reaction rates. In addition, presented results show that DFT computations can be used for reliable prediction of relative reaction rates of vicarious nucleophilic substitution (VNS)-type of reactions. This work provides a valuable tool in the form of structure–thermodynamic property correlation data that can be used in the synthesis design process or QSAR type of analysis, where explicit trends are often required.

---

Introduction

Introduction of substituents into aromatic rings is one of the most important processes in organic chemistry. The most general and widely used process is direct electrophilic substitution in arenes, such as nitration, sulfonation, and Friedel–Crafts reaction.^1–3 These reactions proceed via addition of an electrophilic agent to aromatic rings at positions occupied by hydrogen, followed by fast spontaneous departure of a proton from cationic adducts.^4,5 Therefore, of great interest is the effect of substituents on the rate and orientation of the substitution. The results of the kinetic studies of nitration as the model process of electrophilic substitution laid the foundation for the formulation of the concept of electronic effects of substituents on these reactions.^1–5 Further refining and quantifying the results is the basis of modern physical organic chemistry and the Hammett equation, inter alia.^5–10 These data can be considered quantitative expressions of nucleophilic activities of arenes and even particular positions of aromatic rings.

However, nucleophilic aromatic substitution has long been considered limited to the substitution of halogens located in the activated ortho and para positions of nitroarenes.^11–17 The mechanism of this reaction has been thoroughly studied. It was established that the reaction proceeds via addition of nucleophiles at positions occupied by halogens to form intermediate σ^X adducts, followed by fast departure of the halogen anions. On the basis of kinetic studies, it was found that the addition is the rate limiting step of the substitution. Recently, it was shown that in some cases, the departure of the chlorine anion from the σ^Cl adducts is so fast that the process of addition–elimination proceeds synchronously.¹⁸

Only about 40 years ago, it was unambiguously established that addition of anionic nucleophiles to halonitroarenes proceeds faster at positions occupied by hydrogen to form σ^H adducts than to those occupied by halogens.^19,20 Since hydride anions, contrary to the halogen anions, are unable to depart spontaneously from σ^H adducts, they most often dissociate, and the slower addition at positions occupied by halogens results in S_NAr. Thus, the nucleophilic substitution of halogen is, in fact, a secondary process. Nevertheless, there are a few ways to convert these σ^H adducts into products of the nucleophilic substitution of hydrogen (S_NArH). Variants in the conversion of σ^H adducts, such as the oxidative nucleophilic substitution of hydrogen (ONSH),^21–23 vicarious nucleophilic substitution (VNS),^24–26 and formation of substituted nitrosoarenes, are presently well-established and versatile synthetic methodologies.^27–29 Only when, owing to the nature of nucleophiles and conditions, the fast further conversion of the σ^H adducts cannot proceed, they dissociate, hence allowing the conventional S_NAr of halogens to take place. This corrected mechanistic picture of nucleophilic aromatic substitution was recently presented based on experimental results and quantum chemical calculations.^30–33

Thus, in principle, nucleophilic and electrophilic aromatic substitutions proceed analogously but according to opposite polarity; thus, they are in Umpolung relation.²⁹ It was therefore of great interest to learn the effect of substituents on the electrophilic activities of the nitroaromatic rings. Numerous studies have reported the effects of substituents on the rate of the S_NAr reaction collected in the monographs of Miller and Terrier; however, the results do not reflect virtual activities for two reasons.^11,17 Although the rates of nucleophilic substitution of halogens are equal to the rates of addition, they strongly depend on the nature of the leaving group and particularly the substitution of halogens is secondary reactions preceded by reversible nucleophilic addition at positions occupied by hydrogen. The rates of substitution of halogens can therefore be affected by the equilibrium of the formation of the σ^H adducts.

Considering that the fast initial process between nucleophiles and nitroarenes is an addition at positions occupied by hydrogen, one can expect that the rates of properly selected S_NArH reactions can reflect the rates of nucleophilic addition and be a measure of the effects of substituents on their electrophilic activities.

Some years ago, we determined the relative rates of nucleophilic addition of a model nucleophile, the carbanion of chloromethyl phenyl sulfone to a series of nitroarenes containing various substituents using a competitive VNS reaction under kinetically controlled conditions in relation to the rate of substitution at position ortho of nitrobenzene. These relative rates are virtual measures of the electronic effects of the substituents.³⁴

In this paper, we present a comparison of the experimentally determined relative rates of nucleophilic addition with the quantum chemically calculated thermodynamic parameters of the addition of this model nucleophile at various positions in a series of ortho-, meta- and para-substituted nitroarenes, including halonitroarenes. These parameters include activation of Gibbs free energy (ΔG^#) as well as Gibbs free energy of the reaction (ΔG_R), in which the σ^H adduct is formed. In the last 20 years, quantum chemical calculations have become a routine tool for organic chemists, helping them better understand the reactions they are studying and making possible calculations with different types of spectra and other physicochemical parameters of complex organic molecules.^{31–33,35–39} Owing to hardware limitations, the only practical computational approach for medium-sized molecules is using various methods based on density functional theory (DFT). Unfortunately, the number of established DFT methods is enormous, and there are no general guidelines for selecting a method that is best suited for a given molecular system or reaction. Therefore, in our work, we decided to perform calculations using a series of the most popular methods to observe how the selection of the given DFT method can affect the resulting molecular geometries and their thermodynamic functions.

Qualitative correlation between experimentally determined relative rates of nucleophilic addition with calculated ΔG^# and ΔG_R is the main subject of this work, but we are also interested in two other questions: (i) to what extent is it possible to predict the ratio between ortho and para (relative to the nitro group) substitution products in the case of the ortho-substituted nitrobenzenes and the ratio between three possible products (two ortho and one para) in the case of meta-substituted nitrobenzene derivatives, and (ii) is it possible to predict, which σ^H adducts are stable under reaction conditions, i.e. the equilibrium between substrates and σ^H adducts is shifted in favor of the latter.

Computational methods

All calculations were performed using the Gaussian 16 Rev. A software package.⁴⁰ Starting geometries were created using the GaussView v.6 program.⁴¹ For molecules or ions with more than one possible conformer, the conformational space was searched manually or automatically using Spartan v. 18 software⁴² to find the most stable conformer. The starting geometries of the transition states (TS) were generated by stretching the C–C bond between the benzene ring and the nucleophilic carbon atom of the chloromethyl phenyl sulfone group in the σ^H adduct. The final geometries of the TS were optimized using opt = (TS,noeigen,calcfc) keywords. It was proved that for the optimized TS geometries, one imaginary frequency exists, corresponding to the stretching of the abovementioned C–C bond. The majority of calculations were performed using PBE0 functional (PBE1PBE in Gaussian), B3LYP and ωB97XD. Other functionals that were tested are M062X, PBE1PBE-D3, B3LYP-D3 and APFD. APFD and ωB97XD functionals already contain dispersion corrections in their definitions. PBE1PBE-D3 and B3LYP-D3 notations mean that corrections for dispersion were added using the “EmpiricalDispersion = GD3” keyword. Geometry optimization and frequency calculations were performed using a 6-31+G(d) basis set, while final energies were calculated at the 6-311+G(2d,p) level. To model the liquid phase reactions, the polarizable continuum (PCM) model in the most advanced SMD version was applied with N,N-dimethylformamide (DMF) as the solvent (keyword: scrf = (smd,solvent = N,N-dimethylformamide)).

Results and discussion

As a model nucleophile for this study, we selected the carbanion of chloromethyl phenyl sulfone used previously for most of the mechanistic studies of S_NArH reactions.³⁴ It was also used as a model for the previously reported calculations of the energy profiles of the reaction of nucleophiles with nitrobenzene, p-fluoronitrobenzene and p-chloronitrobenzene.^31–33,43 The present study embraces a large series of nitrobenzenes substituted at various positions ortho, meta and para with various substituents. Calculations were performed in the gas phase and in N,N-dimethylformamide (DMF). In our previous paper, we showed that calculations performed for the gas phase and the DMF solution using the polarizable continuum model (PCM) gave similar relative results concerning the relative values of ΔH and ΔG although absolute values of these functions calculated using the PCM model are significantly lower.³²

An important question was the selection of the DFT method(s) to be used in our study. There is general consent that it is not possible to determine which DFT functional is “the best” in most applications. The most widely used is the B3LYP functional, but in many cases, its performance is questionable. In our previous studies,^31–33,43 we used the PBE0 functional (PBE1PBE in Gaussian) with a 6-31+G(d) basis set for geometry optimization and frequency calculations with 6-311+G(2d,p) basis set for single point energy calculations. In this work, we decided to test a wide range of DFT methods in both the gas phase and the DMF.

para-Substituted nitrobenzenes

The first series of calculations were performed for para-substituted nitrobenzenes because this is the simplest model: σ^H adduct can be formed solely in position ortho to the nitro group. However, the situation is a bit more complex because in the case of the addition of the anion of chloromethyl phenyl sulfone, two diastereoisomeric σ^H adducts can be formed (Scheme 1).⁴⁴

Scheme 1 Formation of diastereoisomeric ortho σ^H adducts of p-substituted nitrobenzenes in reaction with the anion of chloromethyl phenyl sulfone.

In the first part of our work, we modelled the geometries of transition states (TS) and σ^H adducts of the anion of chloromethyl phenyl sulfone to p-chloronitrobenzene in position ortho using our preferred DFT method: PBE1PBE/6-311+G(2d,p)//PBE1PBE/6-31+G(d). p-Chloronitrobenzene was selected as a representative model of p-substituted nitrobenzene derivatives. The results of these calculations, in which full conformational analysis was performed, are presented in Fig. 1. A comparison of the activation Gibbs free energy (ΔG^#) values revealed that the TS (a) leading to adduct (b) with configuration S,R of the stereogenic centers is about 1.1 kcal mol⁻¹ lower than that of TS (c). Additionally, adduct (b) is 1.7 kcal mol⁻¹ more stable than adduct (d). These results agree very well with the established mechanism of the VNS reaction.^{20,27–30,32} In the second step of this reaction, the E2 elimination of HCl induced by a strong base, present in the reaction mixture, occurs, leading to the anion of the final hydrogen substitution product. Such elimination requires the antiperiplanar conformation of the H and Cl atoms. This requirement is fulfilled in the S,R adduct (b) because in its most stable conformation, the H–C–C–Cl dihedral angle is about 177°, a value very close to the optimal 180° required for E2 elimination. In contrast, placing H and Cl atoms in the antiperiplanar conformation in adduct S,S (d) requires several rotations around the C–C and C–S bonds, which makes E2 elimination more difficult.

Fig. 1 Optimized gas-phase geometries of diastereoisomeric transition states (a) and (c) and the resulting σ^H adducts (b) and (d) of the anion of chloromethyl phenyl sulfone to p-chloronitrobenzene in the ortho position. Structure (e) represents the best geometry of the S,S adduct in DMF.

To make sure that the presence of the solvent does not significantly affect the relative stabilities of diastereoisomeric adducts and transition states, the same calculations were performed in N,N-dimethylformamide (DMF) using the PCM-SMD approximation. Again, the S,R and TS adducts were more stable by 0.5 and 2.1 kcal mol⁻¹, respectively. The only important change compared to the gas-phase calculations was the different geometries of the most stable S,S adducts (Fig. 1e). The geometries of both the TS and S,R adducts were very similar in both phases. Owing to these results in all the following calculations, only more stable S,R diastereoisomers and corresponding transition states were considered.

Correlations between the logarithm of the experimental relative reaction rates and the computed thermodynamic functions

In the second stage of our work, the enthalpies and Gibbs free energies of TS of the addition of the anion of chloromethyl phenyl sulfone to a series of p-substituted nitrobenzenes, as well as the energies of the resulting σ^H adducts are calculated. This is made possible to calculate activation enthalpies (ΔH^#) and Gibbs free energies (ΔG^#) as well as reaction enthalpies (ΔH_R) and free energies (ΔG_R) and compare them with the logarithms of the experimental relative rates of reaction of this nucleophile with p-substituted nitrobenzenes, as published by Błażej and Mąkosza.³⁴ All calculations were performed in both the gas phase and DMF using the PCM-SMD solvation model. Because the number of available functionals is enormous and there are no general rules that could help to select an optimal method for the given type of calculations, we decided to test seven popular DFT methods: PBE0 (PBE1PBE in Gaussian), B3LYP, M062X, PBE1PBE-D3, B3LYP-D3, APFD and ωB97XD. The last four functionals include corrections for dispersion, which can be important for calculating transition state geometries and energies. APFD and ωB97XD (Grimme's D2) functionals already contain dispersion correction in their definition, while in the case of the PBE1PBE and B3LYP functionals, corrections for dispersion were added using the “EmpiricalDispersion = GD3” keyword. It must be mentioned that for p-trifluoromethylnitrobenzene, we were not able to find the real energy minimum when calculations were performed in DMF using the PBE1PBE/6-31+G(d) method. Geometry optimization was terminated correctly, i.e. energy minimum was found, but frequency calculations for this geometry showed one imaginary frequency corresponding to the rotation of the CF₃ group. The problem was solved by optimization and frequency calculations using an extended basis set, i.e. 6-311+G(2d,p). The results of the calculations using the PBE0 method are shown in Table 1. Correlations between computed ΔG^# and ΔG_R and the logarithms of experimental relative rate values are presented in Fig. 2.

Table 1 Comparison of the calculated activation Gibbs free energies (ΔG^#) and reaction free energies (ΔG_R) in the gas phase and DMF with experimental relative reaction rates for the reaction of p-substituted nitrobenzenes with the anion of chloromethyl phenyl sulfone using PBE1PBE/6-311+G(2d,p)//PBE1PBE/6-31+G(d) method

–X	log(k) exp.	Gas phase		DMF
		ΔG^#	ΔGR	ΔG^#	ΔGR
a For p-CF3 nitrobenzene in DMF geometry; optimization and frequency calculations were performed using a 6-311+G(2d,p) basis set; see text for explanation.
–iPr	−0.54	8.75	−2.48	21.38	5.91
–t-Bu	−0.44	8.66	−2.61	21.89	6.11
–OCH3	−0.05	7.87	−3.39	21.12	5.87
–H	0.00	8.50	−2.76	21.43	5.68
–OPh	0.43	7.13	−5.41	20.87	5.13
–SMe	0.58	6.45	−5.08	20.60	5.07
–SPh	1.08	5.66	−7.16	20.07	3.65
–COO-iPr	1.20	6.34	−5.70	18.72	3.23
–F	1.70	4.82	−7.22	19.10	2.86
–Cl	2.11	3.95	−8.40	18.96	2.71
–Br	2.18	3.33	−9.54	17.90	1.62
–CF3	2.81	1.99	−10.71	18.41^a	1.45^a
–CN	3.02	−1.12	−14.17	17.16	0.18

---

Fig. 2 Correlations between ΔG^# and ΔG_R and the logarithms of experimental relative rate values for the reaction of p-substituted nitrobenzenes and anion of chloromethyl phenyl sulfone. Upper traces = reaction in DMF; lower traces = reaction in the gas phase. For the reaction in the gas phase, the green correlation line corresponds to the data set without a point for the reaction of p-cyanonitrobenzene (blue dot).

The results presented in Fig. 2 show good correlations between calculated ΔG^# and ΔG_R and logarithms of experimental relative rate values both in the gas phase and in solution in DMF. There is, however, one important discrepancy that can be observed in graphs showing the results of the calculations in the gas phase. The point corresponding to the reaction of p-cyanonitrobenzene is located far below the correlation line, indicating that the computed reaction rate is much higher than the experimental one. Removing it from the correlation increases the R² from 0.92 to 0.97 for ΔG^# and from 0.93 to 0.96 for ΔG_R. This effect is not observed for the calculations performed in DMF. One of the possible rationalizations of this effect considers that in the gas phase, the calculated ΔG^# value for the reaction with p-cyanonitrobenzene is negative, indicating a lack of the activation barrier. The reaction rate in such cases is controlled by diffusion rather than activation energy. In solution, solvation effects increase the ΔG^# value significantly, so all of them are positive and much higher than in the gas phase. It is noteworthy that the entire study is based on the assumption that the second step of the VNS reaction, i.e. E2 elimination is much faster than the addition of nucleophile to the nitroaromatic ring. In the case when addition is fast, the total reaction rate can be controlled by the rate of elimination.

The lack of significant differences for the calculations in the gas phase and in DMF (except for the reaction with p-cyanonitrobenzene) indicates that in this solvent, the anion of chloromethyl phenyl sulfone is relatively free, i.e. it does not form a tight pair with the potassium counterion. This conclusion agrees with experimental studies concerning VNS reactions.^20,27–29 Considering the discussion presented above, it can be concluded that for the prediction of the relative reaction rates of the studied reactions, calculations incorporating the solvation model are recommended. However, for the majority of substrates, gas-phase calculations give comparable results.

The second important question concerns the applicability of the calculated ΔG^# and ΔG_R values for the prediction of the relative reaction rates. For the gas-phase reaction (excluding the result for p-cyanonitrobenzene), the correlation coefficients (R²) are almost the same for correlations with ΔG^# and ΔG_R: 0.97 and 0.96, respectively. In the case of the calculation performed in solution, a significantly better correlation is observed for the free energy of reaction rather than the free energy of activation (R² = 0.96 vs. 0.91). Similar results obtained for ΔG^# and ΔG_R can be rationalized by comparing the structures of the transition states and reaction products. It can be observed in Fig. 1 that the structures of the transition states and σ^H adducts with the S,R configuration, for which calculations were performed, are very similar. This corresponds to the so-called late or product-like transition state. For such reactions, their rates are more dependent on the structures of the products than on the substrates. This result is very useful because it shows that in the case of the reactions of substituted nitrobenzenes with nucleophiles, in which relatively stable σ-adducts are formed, the relative reaction rate can be estimated by modelling the structure of the σ-adduct rather than the transition state, which is a much easier task. A significantly better correlation observed for the ΔG_R rather than ΔG^# in the calculations performed in solution requires further discussion. In our opinion, the most likely rationalization is that quantum chemical calculations can more accurately model transition states in the gas phase than in the solution. In the case of reaction products, the computation accuracy is comparable. Because, as discussed above, the geometries of transition states and corresponding σ-adducts are similar, the final calculations performed in the solution for σ-adducts give more accurate results.

As part of our research, we also attempted to correlate the logarithm of the relative reaction rates not only with Gibbs free energies of activation and reaction but also with respective enthalpy values: ΔH^# and ΔH_R. According to the basic rules of thermodynamics, ΔG^# and ΔG_R should be used to estimate the reaction rate constant and equilibrium constant, respectively. However, computation of Gibbs free energy requires computation of entropy, which can introduce relatively high errors. Provided that the entropy changes in the studied reactions are similar, the correlation with the enthalpy values should be as good as that with Gibbs free energy. The results of this correlation are shown in Fig. S1 in the ESI.† Both correlations for the solution-phase reactions are better than those with respective Gibbs free energy values, while correlations for the gas-phase reactions are lower. These results can be rationalized considering that the computation of entropy is more accurate in the gas phase than in the liquid phase. Moreover, it is not possible to generalize these results, but at least for estimation of the relative rates of addition of the anion of chloromethyl phenyl sulfone to p-substituted nitrobenzenes, the best results are obtained for ΔH_R calculated in DMF solution.

All the results discussed above were obtained using the PBE1PBE/6-311+G(2d,p)//PBE1PBE/6-31+G(d) method, in both the gas phase and the DMF solution. Unfortunately, the selection of the DFT method for solving a given problem is still arbitrary. In our previous studies,^31–33,43 we used the PBE0 functional (PBE1PBE in Gaussian), so this was also our first choice in this work. However, in many recent studies, many other functionals were used. Therefore, we decided to perform our calculations using several of the most popular DFT methods. In Table 2, correlation coefficients R² for the correlation of the logarithm of relative reaction rates with ΔG^# and ΔG_R in the gas phase and with ΔG^#, ΔG_R, ΔH^# and ΔH_R in the DMF solution are collected (full results are presented in ESI†).

Table 2 R ² coefficients for correlation of the logarithm of relative reaction rates of the anion of chloromethyl phenyl sulfone with p-substituted nitrobenzenes with ΔG^# and ΔG_R in the gas phase and with ΔG^#, ΔG_R, ΔH^# and ΔH_R in DMF solution. The best correlations in the gas phase and solution are marked in bold. For information about basis sets and other computational details, see the ‘Computational methods’ section

DFT method	Gas phase		DMF solution
	ΔG^#	ΔGR	ΔG^#	ΔGR	ΔH^#	ΔHR
PBE1PBE	0.92	0.93	0.91	0.96	0.95	0.98
M062X	0.92	0.90	0.91	0.95	0.96	0.98
APFD	0.95	0.91	0.85	0.93	0.91	0.96
PBE1PBE-D3	0.88	0.92	0.94	0.93	0.94	0.97
ωB97XD	0.91	0.90	0.90	0.94	0.95	0.98
ωB97XD (−40 °C)			0.92	0.95	0.95	0.98
B3LYP	0.93	0.95	0.94	0.95	0.95	0.97
B3LYP (−40 °C)			0.95	0.96	0.95	0.97
B3LYP-D3	0.90	0.93	0.90	0.93	0.94	0.97

---

A general conclusion that can be drawn from these results is that there are no significant differences between the majority of the tested DFT methods concerning the correlation between the logarithm of the relative reaction rates of the studied reactions and computed thermodynamic functions in DMF solution. This is true, especially for the correlation with the enthalpy of activation, as well as the enthalpy of reaction. Rather surprisingly, it looks like, for the studied reaction, the best thermodynamic quantity to correlate with the logarithm of the relative reaction rates is the enthalpy of the reaction. For all but one DFT methods, R² values are between 0.97 and 0.98. Correlations with ΔH^# were slightly lower: 0.94 to 0.96 for all DFT methods used. As discussed above, there are two reasons for this effect: (i) much better computational modelling of the reaction product than the transition state and (ii) a late (product-like) transition state for this reaction.

Concerning the better correlation of enthalpy rather than Gibbs free energy values, which was observed for calculations using the PBE0 method, all other DFT methods followed this path. As described above, the most likely rationalization for such results is the relatively low accuracy of the computation of entropy.

The last parameter that can affect the correlation between the logarithm of the relative reaction rates of the studied reactions and the computed thermodynamic functions is temperature. In the original experimental study, all reactions were performed at −40 °C, while calculations were performed at 25 °C as default. There is also experimental evidence that substituted nitrobenzenes can form stable σ^H adducts with the anion of chloromethyl phenyl sulfone at low temperatures.^44,45 To test the temperature effect, we performed calculations of the thermodynamic parameters of the studied reactions using ωB97XD and B3LYP functionals at −40 °C. The results shown in Table 2 indicate that temperature has only a minor effect on the correlation coefficients, so all further calculations were performed under standard conditions.

Reaction profiles

The results of the calculations described in the previous section can also be used to plot the energy profiles of the studied reactions. Such plots showing changes in the Gibbs free energy during the addition of the anion of chloromethyl phenyl sulfone to p-substituted nitrobenzenes calculated using PBE1PBE/6-311+G(2d,p)//PBE1PBE/6-31+G(d) method are presented in Fig. 3.

Fig. 3 Reaction profiles for addition of the anion of chloromethyl phenyl sulfone to p-substituted nitrobenzenes in the gas phase and in solution in DMF obtained using the PBE1PBE/6-311+G(2d,p)//PBE1PBE/6-31+G(d) method.

These results show that in the gas phase, all studied p-substituted nitrobenzenes should form stable σ^H adducts. Relative stability is increased, as expected, with the growing electronegativity of the substituent. Additionally, activation free energies ΔG^# are low, indicating that addition should proceed relatively fast. A different situation is observed in the DMF solution. Activation free energies of addition are significantly higher, and all σ^H adducts are thermodynamically unstable but only by a few kcal mol⁻¹. These results are inconsistent with the experimental observations. It is known that the concentration of σ^H adducts of the nucleophiles to mononitrobenzenes is very low under standard reaction conditions, but such a low concentration is sufficient for very fast E2 elimination. As mentioned above, it is possible to observe stable σ^H adducts to mononitroarenes but only at low temperatures under specific conditions.^44,45

Interesting but rather upsetting conclusions were obtained by comparing the energy profiles of the reactions of selected p-substituted nitroarenes with the anion of chloromethyl phenyl sulfone computed using different DFT methods. Representative examples of three reactions with nitrobenzene, p-chloro- and p-cyanonitrobenzene in DMF solution are presented in Fig. 4.

Fig. 4 Reaction profiles for addition of the anion of chloromethyl phenyl sulfone to nitrobenzene, p-chloro- and p-cyanonitrobenzene in solution in DMF obtained using different DFT methods.

It can be observed that the difference between the results obtained by B3LYP and APFD functionals is about 11 kcal mol⁻¹ for ΔG^# and about 5 kcal mol⁻¹ for ΔG_R. All other functionals we tested provided results between these two extremes. This result clearly shows that at the present stage of the development of DFT methods, it is not possible to reliably estimate the absolute thermodynamic parameters of this type of reaction. Summarizing this part of our work, we can say that it is possible to estimate the relative reaction rate of the VNS reaction by modelling its transition state or the structure of the σ^H adduct. As expected, a better correlation was obtained for the results of the calculations performed in the solution. Surprisingly, better correlations were obtained for ΔH^# and ΔH_R rather than respective Gibbs free energy values, which can indicate problems with accurate calculation of entropy and pure reproduction of subtle solvent-molecule interactions using the implicit solvent model. Finally, it is not possible to estimate, even with low accuracy, the absolute thermodynamic parameters of the addition of the anion of chloromethyl phenyl sulfone to p-substituted nitrobenzenes using DFT methods.

ortho-Substituted nitrobenzenes

The second group of compounds we studied was ortho-substituted nitrobenzenes. In the reaction with the anion of chloromethyl phenyl sulfone, they can form two types of σ^H adducts: in position 4 – para to the nitro group and in position 6 – ortho to the nitro group (Scheme 2).

Scheme 2 Formation of ortho and para σ^H adducts of ortho-substituted nitrobenzenes in reaction with the anion of chloromethyl phenyl sulfone.

In both cases, two diastereoisomeric adducts can be formed. In the case of addition at position 6 (ortho to the nitro group) for all the studied compounds, the S,R diastereoisomer for both TS and the adduct is more stable, similar to the reactions of p-substituted nitrobenzenes (see above). This is true for the reaction in the gas phase as well as for the reaction in DMF. A much more complex situation is observed for the reaction in position 4, i.e. para to the nitro group. An exhaustive conformational search showed that for all studied reactions, the lowest Gibbs free energy has two diastereoisomeric transition states shown in Fig. 5a and c, taking the reaction of 2-chloronitrobenzene as the representative example. For all reactions in DMF and most reactions in the gas phase, geometry (a) is more favorable, but the difference between geometries (a) and (b) is always lower than 0.9 kcal mol⁻¹, which is much lower than the expected accuracy of DFT calculations (for complete results, see ESI†). Consequently, adducts (b) and (d) can be formed with practically equal probability.

Fig. 5 Optimized geometries of the diastereoisomeric transition states (a) and (c) and the resulting σ^H adducts (b) and (d) of the anion of chloromethyl phenyl sulfone to 2-chloronitrobenzene at position 4 (para to the nitro group) computed for the DMF solution. Structure (e) shows the most stable geometry computed for the majority of adducts in the gas phase.

The optimization of the geometries of σ^H adducts in position 4 gave different results for the gas phase and the DMF solution. In DMF, the R,R geometry ((d) in Fig. 5) is preferred for all adducts except for the adduct to 2-nitrobenzoic acid isopropyl ester. The energies of both diastereoisomers are, however, quite similar. The largest difference is 1.1 kcal mol⁻¹, which, considering the accuracy of the DFT calculations, means that both isomers can exist in comparable quantities. For both isomers, the H–C–C–Cl dihedral angle is close to 180°, facilitating the E2 elimination of HCl in the second step of the VNS reaction. In the gas phase, the geometry of most adducts shown in Fig. 5e) is preferred, but the energy differences between the diastereoisomers and conformers are very small. An examination of the geometries of the adducts reveals that in the gas phase, the H–C–C–Cl dihedral angle in the conformers of the lowest energies differs significantly from the optimal 180°. In conclusion, in the case of the addition in position 4, calculations performed considering the presence of the solvent give more reliable results.

The results of the calculations performed in the gas phase and the DMF solution using PBE1PBE/6-311+G(2d,p)//PBE1PBE/6-31+G(d) method are presented in Table 3 and Fig. 6. Plots in Fig. 6 are drawn separately for addition in positions 4 and 6. Similarly, to the results obtained for p-substituted nitrobenzenes, a good correlation between calculated ΔG^# and ΔG_R and logarithms of experimental relative rate values in both the gas phase and the DMF solution is observed. There is, however, an important difference between the correlations of ΔG^# and ΔG_R. In the case of the correlation of the reaction between Gibbs free energy and the log(k_exp), the plots for the reactions in positions 4 and 6 practically overlap. Such a situation should be expected provided that calculations reproduce an addition in both positions with the same accuracy. In the case of activation, the Gibbs free energy plots for the reactions in positions 4 and 6 still show good linearity, but the correlation lines, although almost parallel, are shifted by about 3 kcal mol⁻¹ for the reactions in the gas phase and by less than 2 kcal mol⁻¹ for the reactions in DMF. These results demonstrate that our calculation method shows different systematic errors when applied to model the transition states of additions in positions ortho and para to the nitro group.

Table 3 Comparison of the calculated activation Gibbs free energies (ΔG^#) and reaction free energies (ΔG_R) in the gas phase and DMF with experimental relative reaction rates for the reaction of o-substituted nitrobenzenes with the anion of chloromethyl phenyl sulfone in positions 4 and 6 using PBE1PBE/6-311+G(2d,p)//PBE1PBE/6-31+G(d) method

–X	Reaction in pos. 4					Reaction in pos. 6
	log(k) exp.	Gas phase		DMF		log(k) exp.	Gas phase		DMF
		ΔG^#	ΔGR	ΔG^#	ΔGR		ΔG^#	ΔGR	ΔG^#	ΔGR
–H	−0.36	5.70	−1.73	20.21	5.62	0.00	8.50	−2.76	21.43	5.68
–COO-iPr	2.27	3.20	−6.42	18.35	3.27	2.71	5.13	−7.75	18.70	2.12
–F	3.18	1.96	−5.95	17.62	1.98	3.04	5.40	−6.38	18.85	2.28
–Cl	2.17	2.33	−6.93	18.57	2.10	2.15	5.43	−6.73	19.26	2.10
–Br	1.82	2.50	−6.87	18.72	4.11	1.93	5.32	−7.83	20.69	1.84
–CF3	4.13	−0.24	−9.58	17.21	0.40	4.79	2.75	−9.90	18.08	0.52
–CN	7.09	−3.22	−12.74	15.19	−1.40	7.44	−0.48	−13.52	17.40	−1.71

---

Fig. 6 Correlations between ΔG^# (left) and ΔG_R (right) and the logarithms of experimental relative rate values for the reaction of o-substituted nitrobenzenes with the anion of chloromethyl phenyl sulfone in positions 4 (hollow circles) and 6 (filled circles) using the PBE1PBE/6-311+G(2d,p)//PBE1PBE/6-31+G(d) method. Upper traces = reaction in DMF; lower traces = reaction in the gas phase.

To check the influence of the selected computational method, we repeated the calculations in the DMF solution using the B3LYP and ωB97XD functionals with the same basis sets, which were used for calculations using the PBE0 method. The results are shown in Fig. S2–S4 in ESI.† Probably, the most important conclusion that can be drawn from them is that all three DFT methods provide a very good or at least reasonable correlation between the experimental and computed data. This is true for ΔG^# and ΔG_R, as well as for respective enthalpy values. The best correlation was obtained for ΔH^# calculated using the B3LYP functional: 0.98 both for addition in positions 4 and 6. Concerning the ΔH_R calculated using this functional, an almost perfect correlation was obtained for addition in position 4 (R² = 0.997) and lower for addition in position 6 (R² = 0.97).

Examination of the plots from Fig. S2–S4 (ESI†) shows that on the majority of them, correlation lines for addition in positions 4 and 6 are parallel but shifted by up to 2 kcal mol⁻¹. Unfortunately, this shift has different signs when comparing different methods and thermodynamic quantities. For example, in the case of ΔG^# and ΔH^#, higher values for addition in position 6 are obtained using PBE0 and B3LYP functionals, while in the case of ωB97XD functional, the opposite is true. The only general but not very useful conclusion that can be drawn is that each DFT method introduces a different systematic error in modelling addition in positions 4 and 6. This also means that at the present stage of the development of the computational methods, it is not possible to predict the ortho to para ratio of the nucleophilic addition reaction to ortho-substituted nitrobenzene derivatives. The subtle factors affecting the kinetic and thermodynamic nature of the process, such as actual solvent–molecule interactions, temperature dependence and intra-atomic interaction between atoms in ortho and meta-oriented atomic systems, are difficult to reproduce by applying DFT methods when combined in one molecule. To achieve this level of accuracy the results of the calculations should be accurate up to about 0.1 kcal mol⁻¹, which is impossible for such molecular systems.

meta-Substituted nitrobenzenes

The third group of compounds studied in this work is meta-substituted nitrobenzene derivatives. For these compounds, the substitution of hydrogen is possible in three different positions (Scheme 3).

Scheme 3 Formation of 2, 4 and 6 σ^H adducts of meta-substituted nitrobenzenes in reaction with the anion of chloromethyl phenyl sulfone.

Conformational analysis performed for the reactions in DMF solution showed that in the case of addition in position 6 (ortho to nitro group and para to substituent X), only the S,R diastereoisomeric adduct predominates for all studied compounds (Fig. 7h), as observed for the reactions of ortho and para-substituted nitrobenzenes. It must be noted, however, that the energy difference between transition states, leading to S,R and S,S diastereoisomeric σ^H adducts, is very low for the majority of studied reactions. This result indicates that the formation of both diastereoisomeric adducts is possible.

Fig. 7 Optimized geometries of the diastereoisomeric transition state and resulting σ^H adducts of the anion of chloromethyl phenyl sulfone to 3-chloronitrobenzene: in position 2 (a)–(d), 4 (e) and (f) and 6 (g) and (h) computed for the DMF solution. Differences between ΔG in kcal mol⁻¹.

A similar situation is observed for addition in position 2 (ortho to both substituents), where the S,R diastereoisomer is most stable for almost all adducts (Fig. 7c). The only exception is the 3-trifluoromethyl derivative; however, in this case, the energy difference between S,R and S,S diastereoisomeric σ^H adducts is negligible (0.03 kcal mol⁻¹). Concerning the addition in position 4 (para to a nitro group and ortho to substituent X), a less sterically hindered S,S isomer predominates. Its stereochemistry also facilitates the E2 elimination of HCl owing to the antiperiplanar position of hydrogen and chlorine atoms (Fig. 7f).

Similar to two other isomers, for meta-substituted nitrobenzenes, calculations were performed in the gas phase and DMF solution using PBE1PBE/6-311+G(2d,p)//PBE1PBE/6-31+G(d) method. The results are presented in Table 4 and Fig. 8. It is noteworthy that for some substrates, only two or even one substitution product was observed (empty fields in Table 3). Plots are drawn separately for addition in positions 2, 4 and 6. A good correlation between calculated ΔG^# and ΔG_R and logarithms of the experimental relative rate values in both the gas phase and DMF solution is observed for addition in positions 4 and 6 and much worse for addition in position 2. This last result is not surprising because this position is the most sterically hindered, so many subtle effects, which are not correctly described by DFT calculations, can affect the final results. A higher correlation for ΔG^# than for ΔG_R supports this supposition because in TS, the distance between reacting species is larger than in adduct, so steric effects are weaker.

Table 4 Comparison of the calculated activation Gibbs free energies (ΔG^#) and reaction free energies (ΔG_R) in the gas phase and DMF with experimental relative reaction rates for the reaction of m-substituted nitrobenzenes with the anion of chloromethyl phenyl sulfone in positions 2, 4 and 6 using the PBE1PBE/6-311+G(2d,p)//PBE1PBE/6-31+G(d) method

–X	Reaction in pos. 2					Reaction in pos. 4					Reaction in pos. 6
	log(k) exp.	Gas phase		DMF		log(k) exp.	Gas phase		DMF		log(k) exp.	Gas phase		DMF
		ΔG^#	ΔGR	ΔG^#	ΔGR		ΔG^#	ΔGR	ΔG^#	ΔGR		ΔG^#	ΔGR	ΔG^#	ΔGR
a No TS could be found for this reaction.
–Me	−1.02	9.70	0.13	22.07	9.28	−0.49	7.50	0.60	21.06	5.58	−1.66	10.32	−0.49	22.07	7.15
–OMe	−0.48	10.22	−1.59	22.33	6.21	0.79	7.23	−1.52	20.52	2.90	−1.56	10.45	−0.33	22.36	6.99
–H	0.00	8.50	−2.76	21.43	5.68	−0.36	5.70	−1.73	20.21	6.94	0.00	8.50	−2.76	21.43	5.68
–COOMe						6.88	1.41	−14.56	17.15	−3.05	8.27	0.96	−15.03	16.51	−2.96
–F	2.89	6.87	−4.66	19.48	4.16	3.26	3.36	−5.28	17.94	1.65	1.19	6.51	−4.22	20.06	5.07
–Cl	6.38	6.43	−5.24	19.54	4.41						4.70	4.07	−7.94	18.67	2.12
–Br	5.91	6.36	−5.46	19.88	4.85						4.47	3.14	−9.12	18.37	1.31
–CF3											8.52	0.05	−15.05	15.74	−2.53
–CN	8.29		−14.86	16.71	−1.92	9.21	−5.02	−19.68	13.89	−6.77	7.17	−4.59	−20.61	14.54	−6.35

---

Fig. 8 Correlations between calculated ΔG^# and ΔG_R values and the logarithms of experimental relative rate values for the reaction of m-substituted nitrobenzenes with the anion of chloromethyl phenyl sulfone in position 2 (circles), 4 (squares) and 6 (triangles). Upper traces = reaction in the DMF; lower traces = reaction in the gas phase. Black correlation lines and corresponding R² values regarding an addition in position 6 without 3-cyanonitrobenzene. Calculations were performed using the PBE1PBE/6-311+G(2d,p)//PBE1PBE/6-31+G(d) method in the DMF.

meta-Substituted nitrobenzene derivatives differ from the two other isomers because all three possible positions in which nucleophilic addition can occur are oriented ortho or para not only with respect to the nitro group but also to the substituent X. This means that electron withdrawing substituents can significantly affect the rate of addition. The best example is the addition of position 6 to m-cyanonitrobenzene. Computed ΔG^# and ΔG_R values in both the gas phase and the DMF solution are evidently too low compared to the experimental reaction rate value. This means that the reaction rate is controlled by terms other than just the activation of Gibbs free energy, most likely by diffusion. After the removal of the points corresponding to this reaction from the plots in Fig. 8, the correlation coefficients grow very significantly (green and black numbers illustrated in Fig. 8). In contrast to the similar effect observed for the addition of the anion of chloromethyl phenyl sulfone to p-cyanonitrobenzene, in the case of m-cyanonitrobenzene, the effect is also observed in the DMF solution. What makes the results more difficult to rationalize is the lack of a significant discrepancy between the computed and experimental results for addition in position 4 in both the gas phase and, especially, in the DMF.

Calculations performed in the DMF solution using the B3LYP and ωB97XD functionals with the same basis sets, which were used for the calculations using the PBE0 method, showed that for meta-substituted nitrobenzene derivatives, all three functionals gave comparable results (see ESI†). The best correlations with all four computed thermodynamic parameters (ΔG^#, ΔG_R, ΔH^# and ΔH_R) are observed for the addition in position 4 and worst for addition in position 2, which is the most sterically crowded. In general, a rather disappointing conclusion is that for relative rates of addition of a nucleophile to meta-substituted nitrobenzenes, computational predictions are of rather low quality and are differentially reproduced by applying each DFT method. This can have its origin in pure electronic description of the complex system with key thermodynamic and kinetic factors hidden behind weak and subtle solvent-molecule and intra-atomic interactions dedicated to structure description; dispersion-corrected functionals can reproduce properly.^46–49 It is assumed that the electron energy of the systems, ergo, the proper molecular geometries in which the meta-driven substitution electron effect is the strongest, is not equally described by each hybrid and dispersion-corrected functionals.

Conclusions

Modelling of the rate-limiting and key, first step of the VNS reaction between ortho-, meta- and para-substituted nitrobenzene derivatives and the anion of chloromethyl phenyl sulfone using various DFT methods led to the following conclusions.

(1) According to a well-established mechanism of the VNS reaction in the second step, the base-induced β-elimination of HX (X = Cl in the case of the reactions studied in this work) from the preliminarily formed σ^H adduct takes place. This requires antiperiplanar conformation of H and X, i.e. the dihedral angle H–C–C–X in σ^H adduct should be close to 180°. In all reactions studied in this work, two diastereoisomeric σ^H adducts can be formed (the only exception is the addition of position 4 to nitrobenzene). Our calculations showed that in the overwhelming majority of reactions, the diastereoisomer, in which the H–C–C–Cl dihedral angle in the conformer with the lowest energy is close to 180°, is more stable. This result gives some rationalization to the observation that the VNS reaction is usually very fast, even for nitroarenes with relatively low electrophilicity.

(2) Computations using different DFT methods show that all studied functionals give good correlations between the computed thermodynamic parameters and the logarithms of the experimental relative reaction rates. These results show that DFT computations can be used for reliable prediction of the relative reaction rates of the VNS reactions of the substituted nitrobenzene derivatives. The electron withdrawing effect of different functional groups located on the aromatic ring has been reproduced correctly in the presented trends. However, this statement is true only for comparing the reaction rates of addition in the same position of the nitrobenzene derivatives bearing the substituent in the same position relative to the nitro group.

(3) Comparison of the logarithms of the experimental relative reaction rates with computed ΔG^#, ΔG_R, ΔH^# and ΔH_R in DMF solutions shows that all four thermodynamic parameters give good correlations, but for most cases, correlation with ΔH_R is the best.

(4) Better results were obtained for calculations performed in DMF solution, so this approach is recommended.

(5) It is possible to predict separately the quantitative and qualitative trends of the relative rates of substitution in different positions in ortho- and meta- and para-substituted nitrobenzene derivatives, in which the intramolecular interactions between atoms are structurally equal within the same reaction center.

(6) It is not possible to compare absolute and relative kinetic parameters between ortho-, meta- and para-substituted nitrobenzenes because differences in subtle intra-molecular interactions between nucleophile and electrophile are not well reproduced by applying tested DFT methods.

(7) It is not possible to reliably estimate the exact (nominal) thermodynamic parameters of this type of reaction because different DFT methods give different results with differences, reaching up to 15 kcal mol⁻¹. All tested DFT methods with and without dispersion correction seem to reproduce molecular geometries and final energies with comparable systematic errors. The lack of experimental data concerning the stability of anionic σ^H-adducts of nucleophiles to mono-substituted nitrobenzenes makes it impossible to determine the best DFT method empirically.

(8) Presented results in the form of structure–activity trends are valuable tools to analyze the electron withdrawing effects induced by different functional groups in the molecular system within one chosen type of ortho-, meta- and para-oriented molecules.

Author contributions

KB, MM, PS, WD contributed equally. All authors read and approved the final manuscript.

Data availability

The data that supports the findings of this study are available in the ESI† of this article.

Conflicts of interest

The authors declare no competing financial interest.

Acknowledgements

This work has been financed by the National Science Centre, Poland, grant SONATA 17 no. 2021/43/D/ST4/01679. We express our thanks to the Wroclaw Center for Networking and Supercomputing (WCSS) and Interdisciplinary Centre for Mathematical and Computational Modeling (ICM) in Warsaw (grant no. G50-2) for providing computer time and facilities. We gratefully acknowledge Polish high-performance computing infrastructure PLGrid (HPC Centers: ACK Cyfronet AGH) for providing computer facilities and support within computational grant no. PLG/2024/017834.

References

1. R. Taylor, Electrophilic Aromatic Substitution, John Wiley & Sons, Inc., New York, 1990.
2. J. March, Advanced organic chemistry: reactions, mechanisms, and structure, McGraw-Hill, New York, 1977.
3. L. Stock, in A Classic Mechanism for Aromatic Nitration, Progress in Physical Organic Chemistry, ed. R. W. Taft, 1976, pp. 21–47.
4. G. Olah, R. Mahotra and S. C. Narang, Nitration: Methods and Mechanism, John Wiley & Sons, Inc., New Jork, 1989.
5. C. Ingold, Structure and Mechanism in Organic Chemistry, Cornell University Press, 1969.
6. L. P. Hammett, J. Am. Chem. Soc., 1937, 59, 96–103.
7. C. Hansch, A. Leo and R. W. Taft, Chem. Rev., 1991, 91, 165–195.
8. J. Mortier, Arene Chemistry: Reaction Mechanisms and Methods for Aromatic Compounds, John Wiley & Sons, 2016.
9. O. Exner and K. Zvára, J. Phys. Org. Chem., 1999, 12, 151–156.
10. C. D. Ritchie and W. F. Sager, Prog. Phys. Org. Chem., 1964, 2, 323–400.
11. J. Miller, Aromatic Nucleophilic Substitution, Elsevier, New York, 1968.
12. S. Caron and A. Ghosh, Practical Synthetic Organic Chemistry, John Wiley & Sons, Inc, 2011, pp. 237–253.
13. E. Buncel, M. R. Crampton, M. J. Strauss and F. Terrier, Electron-Deficient Aromatic- and Heteroaromatic-Base Interactions, Elsevier, Amsterdam, 1984.
14. J. F. Bunnett, Q. Rev., Chem. Soc., 1958, 12, 1–16.
15. J. F. Bunnett and J. J. Randall, J. Am. Chem. Soc., 1958, 80, 6020–6030.
16. J. F. Bunnett and R. E. Zahler, Chem. Rev., 1951, 49, 273–412.
17. F. Terrier, Nucleophilic Aromatic Displacement, Verlag Chemie, Weinheim, 1991.
18. A. J. J. Lennox, Angew. Chem., Int. Ed., 2018, 57, 14686–14688.
19. J. Goliński and M. Makosza, Tetrahedron Lett., 1978, 19, 3495–3498.
20. M. Mąkosza and J. Winiarski, Acc. Chem. Res., 1987, 20, 282–289.
21. M. Mąkosza and K. Staliński, Chem. - Eur. J., 1997, 3, 2025–2031.
22. M. Makosza, K. Stalinski and C. Klepka, Chem. Commun., 1996, 837–838.
23. D. Sulikowski and M. Mąkosza, Eur. J. Org. Chem., 2010, 4218–4226.
24. J. Goliński and M. Makosza, Synthesis, 1978, 823–825.
25. B. Mudryk and M. Makosza, Tetrahedron, 1988, 44, 209–213.
26. M. Makosza, Synthesis, 1991, 103–111.
27. M. Makosza and T. Glinka, J. Org. Chem., 1983, 48, 3860–3861.
28. M. Mąkosza and K. Wojciechowski, Chem. Rev., 2004, 104, 2631–2666.
29. M. Mąkosza, Chem. - Eur. J., 2020, 26, 15346–15353.
30. M. Mąkosza, Synthesis, 2017, 3247–3254.
31. K. Błaziak, W. Danikiewicz and M. Mąkosza, Molecules, 2020, 25, 4819.
32. K. Błaziak, W. Danikiewicz and M. Mąkosza, J. Am. Chem. Soc., 2016, 138, 7276–7281.
33. K. Błaziak, M. Mąkosza and W. Danikiewicz, Chem. - Eur. J., 2015, 21, 6048–6051.
34. S. Błażej and M. Mąkosza, Chem. - Eur. J., 2008, 14, 11113–11122.
35. T. M. El-Gogary and A. M. El-Nahas, THEOCHEM, 2008, 851, 54–62.
36. M. Govindarajan, K. Ganasan, S. Periandy and S. Mohan, Spectrochim. Acta, Part B, 2010, 76, 12–21.
37. Q.-L. Xu, H. Gao, M. Yousufuddin, D. H. Ess and L. Kürti, J. Am. Chem. Soc., 2013, 135, 14048–14051.
38. D. Avcı, S. Bahçeli, Ö. Tamer and Y. Atalay, Can. J. Chem., 2015, 93, 1147–1156.
39. D. Antoniak, B. Pałuba, T. Basak, K. Błaziak and M. Barbasiewicz, Chem. - Eur. J., 2022, 28, e202201153.
40. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg Williams, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian, Wallingford, CT, 2016.
41. R. Dennington, T. Keith and J. Millam, GaussView ver. 5, Semichem Inc., Shawnee Mission, KS, 2009.
42. L. F. M. Y. Shao, Y. Jung, J. Kussmann, C. Ochsenfeld and S. T. Brown, et al. , Phys. Chem. Chem. Phys., 2006, 8, 3172.
43. K. Błaziak, P. Sendys and W. Danikiewicz, ChemPhysChem, 2016, 17, 850–858.
44. T. Lemek, M. Mąkosza, D. S. Stephenson and H. Mayr, Angew. Chem., Int. Ed., 2003, 42, 2793–2795.
45. B. Chiavarino, M. E. Crestoni, S. Fornarini, F. Lanucara, J. Lemaire and P. Maître, Angew. Chem., Int. Ed., 2007, 46, 1995–1998.
46. M. Wang, X. He, M. Taylor, W. Lorpaiboon, H. Mun and J. Ho, J. Chem. Theory Comput., 2023, 19, 5036–5046.
47. E. Paulechka and A. Kazakov, J. Chem. Theory Comput., 2018, 14, 5920–5932.
48. A. S. Flemming, B. C. Dutmer and T. M. Gilbert, ACS Omega, 2023, 8, 14160–14170.
49. V. Ojea and M. Ruiz, Dalton Trans., 2024, 53, 8662–8679.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp00464g

---