Insight from energy decomposition analysis on a hydrogen-bond-mediated mechanism for on-water catalysis

Abstract

Many classes of organic reactions exhibit a remarkable increase in reaction rates when they occur at the water–organic interface. Although this observed ‘on-water’ catalysis has been extensively studied, the suggested mechanisms still do not explain some of the experimental findings. The mechanism proposed by Jung and Marcus (JACS 129, 5492 (2007)) involves stabilising the transition-state (TS) complex via H-bonds to ‘dangling’ interfacial water molecules. Although the reactants also experience H-bonding to interfacial water molecules in the reactant configuration, it has been argued that the H-bonds are enhanced, in terms of number and strength, in the TS. Therefore, the observed decrease in activation energy has been attributed to this preferential enhancement of H-bonds which leads to a more pronounced TS stabilisation. We employ energy decomposition analysis using the method of absolutely localised molecular orbitals to study this proposition. We find that H-bonds to interfacial water molecules are equivalent in the TS and reactant configurations. Nevertheless, these H-bonds result in significantly enhanced charge-transfer between the reactants in the TS complex, which rationalises the decrease in activation energy.

GRAPHICAL ABSTRACT

Keywords:

• ALMO
• energy decomposition analysis
• on-water catalysis

View correction statement:

Correction

1. Introduction

The field of on-water catalysis added a greener and economic dimension to the current toolbox of organic synthesis. The term ‘on-water’ refers to reactions in which the transition-state (TS) occurs at the organic–water interface on the organic side [1]. This is to differentiate in-water catalysis, which had been previously investigated and refers to reactions occurring in the bulk of the solution [2]. Over the past two decades, many experimental and theoretical studies contributed to the current understanding of the on-water phenomenon [2]. There are three main observations associated with on-water catalysis [3]. Firstly, heterogeneity of the reaction mixture is essential. Secondly, the interface must be with water. Thirdly, some reactions tend to show a large isotope effect, while others show a very mild one. In-water catalysis has been mainly attributed to hydrophobic effects and the cohesive energy density of water [2,4]. On the other hand, on-water catalysis has been rationalised using two mechanisms. The first mechanism involves lowering of the activation energy via H-bonds from dangling OH groups at the water interface [4]. The presence of these dangling OH groups with extra liability to H-bonding has been supported by various experimental and theoretical studies [2,4]. The second catalytic mechanism, proposed by Beattie and coworkers, involves proton-transfer to the organic reactants [3].

The cycloaddition reaction between quadricyclane (1) and dialkyl azodicarboxylates is a good example to investigate the mechanism of on-water catalysis (Figure 1). The cycloadditions involving 1 with dimethyl azodicarboxylate (DMAzD, 2a) or diethyl azodicarboxylate (DEAzD, 2b) show different catalytic rates and isotope effects [2]. Using a TS model for the reaction of 1 and 2a that forms H-bonds with the dangling OH groups at the water interface, Jung and Marcus calculated a drop in activation energy of 7.5 kcal/mol. This value is very close to the experimental value (5.4 kcal/mol) for the on-water cycloaddition of 1 and 2b [5]. Nevertheless, catalysis via H-bonding does not explain the strong isotope effect of the reaction of 1 and 2a [3]. Domingo et al. computed a similar fall in activation energy (8.5 kcal/mol) for the cycloaddition of 1 and 2a albeit adopting a polar mechanism [6]. The mechanism by Domingo et al. is more compatible with the proton-transfer explanation of Beattie and coworkers [2]. The proton-transfer mechanism explains the isotope effect observed in many on-water catalysed reactions and is consistent with the fact that most of the reactions that are catalysed on-water are also subject to acid-base catalysis [3]. Being a hot topic of research, there are a plethora of experimental and theoretical studies endorsing either H-bond or proton-transfer explanations of on-water catalysis [2].

Figure 1. Cycloaddition reaction of quadricyclane (1) and dialkyl azodicarboxylates (2) and three derivatives (2a, 2b and 2c).

In particular, Jung and Marcus concluded that the on-water reaction between 1 and 2a was catalysed due to extra stabilisation of the TS. In their reaction model, 2a makes three H-bonds with water in the TS configuration, while it makes only two H-bonds in the reactant one. They also computed partial charges on the nitrogens and carbonyl oxygens of 2a in the TS and reactant configurations. They deduced that 2a also forms stronger H-bonds to interfacial water molecules in the TS configuration which had more electronegative atoms than the reactant counterpart. Accordingly, they predicted that replacing the nitrogens in DMAzD (2a) with an olefin (2c) or an acetylene group with lower capability of H-bonding would result in inefficient catalysis [4]. Worth noting, Beattie et al. showed experimentally that the cycloaddition of cyclopentadiene and 2c actually exhibits on-water catalysis [3]. Later, Karhan et al. studied the same cycloaddition of 1 and 2a using second-generation Car–Parrinello molecular dynamics [7]. Only a marginal increase of H-bonds has been noticed at the TS. The focus of this work is thus to investigate the relative strengths of H-bonds between 2a and interfacial water molecules in the reactant and TS models.

2. Computational methods

We performed energy decomposition analysis (EDA) of the reactants and TS models for the reaction of 1 and 2a. The method of absolutely localised molecular orbitals (ALMO), available in the CP2K suite (version 7.0)[8], was used to decompose the interaction energy into frozen density, polarisation and charge-transfer (CT) components. Recent reviews are already available on EDA [9,10]. We briefly outline the ALMO EDA method, while the interested reader can find more details and applications in the cited literature [11–14].

The ALMOs on a given fragment, x, are used to compute the state ${\mathrm{\Psi }}_x$ which is infinitely separated from all other fragments. The occupied ALMOs of each fragment are then used to yield an antisymmetrised wave functions ${\mathrm{\Psi }}_0$ of the complex composed of all fragments. At that stage, the ALMOs on a given fragment are not relaxed, except by the changes needed to satisfy the Pauli exclusion principle. The frozen energy contribution is defined as the difference between the SCF energies of those two states: (1) $\mathrm{\Delta }{E}_{\mathrm{FRZ}}=E({\mathrm{\Psi }}_0)-\sum _{x}E({\mathrm{\Psi }}_x)$(1) which means that it includes the Coulomb term as well as an exchange-correlation term. In the next step, the ALMOs on a given fragment are allowed to relax in the field of all other fragments to yield the state ${\mathrm{\Psi }}_{\mathrm{POL}}$ The ALMOs of ${\mathrm{\Psi }}_{\mathrm{POL}}$ are still localised at their fragments. The polarisation energy is then defined as: (2) $\mathrm{\Delta }{E}_{\mathrm{POL}}=E({\mathrm{\Psi }}_{\mathrm{POL}})-E({\mathrm{\Psi }}_0)$(2) The final portion of interaction energy is defined as the difference between the energy of the fully delocalised state ${\mathrm{\Psi }}_{\mathrm{DEL}}$ and the polarised state: (3) $\mathrm{\Delta }{E}_{\mathrm{CT}}=E({\mathrm{\Psi }}_{\mathrm{DEL}})-E({\mathrm{\Psi }}_{\mathrm{POL}})$(3) where $\mathrm{\Delta }{E}_{\mathrm{CT}}$ is the energy lowering due to the total electron transfer from the occupied ALMOs on one fragment to the virtual orbitals on the others (forward donation), in addition to the back donation due to repolarisation.

Our ALMO calculations were conducted at the M06-2x/aug-cc-pVTZ level of theory using the Gaussian-augmented plane wave all-electron approach [8,15–17]. Five setups where initially used in total, as shown in Figure 2: 2a with 3 water molecules (the 3 water molecules are collectively referred to as ‘ 3’) in the reactants configuration (i), and the TS comprised of the coordinates of 1, 2a and 3(ii). The coordinates for i and ii were taken from the study by Jung and Marcus [4] and are provided in Table S1 along with key geometric elements in Figure S1 in the electronic supplementary information (ESI). The three remaining setups were derived from ii to get insight on the EDA of the interaction energy between a pair of 1, 2a and 3 in the absence of the third. This way, the role of including water molecules in the transition state complex can be discerned. The coordinates of 1 were deleted to investigate the interaction between 2a and 3 in the TS conformation (iii). Similarly, the coordinates of 3 were deleted from ii to study the interaction between 1 and 2a maintaining the TS conformation (iv). Finally, the coordinates of 2a were deleted from (ii) to investigate possible interactions between 1 and 3 in the TS conformation (v). Alternative partitioning of the system, taking the three water molecules as a single fragment and thus excluding the inter-water energy components, is provided in Figure 3.

Figure 2. The five systems studied by EDA: (i) reactant conformation of 2a and 3, (ii) TS, (iii) TS without 1, (iv) TS without 3 and (v) TS without 2a. The CT between a given pair of fragments is indicated by an arrow and the corresponding value is given in kJ/mol. Water molecules are collectively referred to as 3, but each water molecule is treated here as a separate fragment.

Figure 3. An alternative partitioning of the studied systems in Figure 2, but combining the water molecules in one fragment: (i.a) reactant conformation of 2a and 3, (ii.a) TS, (iii.a) TS without 1 and (v) TS without 2a. The CT between a given pair of fragments is indicated by an arrow and the corresponding value is given in kJ/mol.

In addition to ALMO, we also did an analysis of configurations i, ii, iii, iv and v using natural bond orbitals (NBO) [18,19]. The NBO computations were done at M06-2x/cc-pVDZ, as implemented in the Gaussian software [16,20]. The relevant NBO results are provided in Figure S2 in the ESI.

3. Results

The five initial configurations are illustrated in Figure 2 with detailed values for CT, while the total results for the EDA of the five systems are provided in Table 1. The CT energies between each pair of 1, 2a and 3 are extracted from Figure 2 and provided in Table 2. In all configurations in Figure 2, each molecule represents a fragment. Hence, the values in Table 2 also includes interaction energy components between water molecules. The EDA results for the alternative partitioning, taking the water molecules as one fragment, are similarly provided in Figure 3 and Tables 3 and 4.

Table 1. EDA results in kJ/mol for the five configurations illustrated in Figure 2.

Energy component	i	ii	iii	iv	v
Frozen	−45.3	78.0	−26.2	102.7	−6.1
Polarisation	−28.4	−70.8	−28.8	−29.6	−13.0
Charge transfer	−22.8	−135.2	−19.0	−84.7	−11.8

ALMO computations were done at the M06-2x/aug-cc-pVTZ level of theory.

Table 2. CT (in kJ/mol) between the fragments in each of the five configurations illustrated in Figure 2.

Involved fragments	i	ii	iii	iv	v
3 and 2a	−12.4	−12.0	−12.0	–	–
3 and 1	–	−5.3	–	–	−4.9
2a and 1	–	−110.6	–	−84.7	–

ALMO computations were done at the M06-2x/aug-cc-pVTZ level of theory.

Table 3. EDA results in kJ/mol for the four configurations illustrated in Figure 3.

Energy component	i.a	ii.a	iii.a	v.a
Frozen density	−33.0	85.3	−18.5	1.2
Polarisation	−17.6	−63.2	−21.4	−5.0
Charge transfer	−12.4	−128.3	−12.0	−4.9

ALMO computations were done at the M06-2x/aug-cc-pVTZ level of theory.

Table 4. CT (in kJ/mol) between the fragments in each of the four configurations illustrated in Figure 3.

Involved fragments	i.a	ii.a	iii.a	v.a
3 and 2a	−12.4	−11.9	−12.0	–
3 and 1	–	−5.4	–	−4.9
2a and 1	–	−110.9	–	–

ALMO computations were done at the M06-2x/aug-cc-pVTZ level of theory.

To establish a reference to compare the magnitudes of energy contributions, we did EDA for the 3 water molecules in the TS configuration at the same level of theory, but isolated from 1 and 2a. The values for frozen, polarisation and CT contributions are −7.6, −7.8 and −6.9 kJ/mol, respectively. As previously reported, frozen, polarisation and CT contributions are of equivalent importance in H-bonds [12].

4. Discussion

As mentioned above, Jung and Marcus concluded that H-bonding is enhanced in the TS based on the expected electrostatics between the computed partial charges in both configurations. Interestingly, the frozen density interaction between 2a and 3 is more enhanced in i than in iii by −15 kJ/mol, see Table 3. Therefore, based on the electrostatic term alone, 2a is expected to form stronger H-bonds to interfacial waters in the reactant conformation than in the TS. This is an example of how conclusions based on population analysis could be misleading. However, water decreased the frozen density repulsion from 103 kJ/mol (between 1 and 2a in iv) to 78 kJ/mol (for the total system, ii). This decrease can be attributed to favourable frozen density attraction within water molecules (−7.6 kJ/mol for the isolated water trimer, see above) and through H-bonds to 2a (−18.5 kJ/mol for iii.a, see Table 3).

We further analyse the remaining two terms that contribute to the H-bond strength. As summarised in Table 2, the CT between 2a and 3 is equivalent in both reactant (i) and TS (ii) configurations. The value is also equivalent when 2a and 3 are isolated from 1 in the TS configuration (iii). Clearly, the CT term does not indicate stronger H-bonding between 2a and 3 in the TS. The same conclusion is reached when the three water molecules are treated as one fragment, see Table 4. We actually notice additional CT between the water molecules and 1 which might contribute, albeit by a tiny fraction, to the energy lowering of the TS. More importantly, the presence of water molecules increases CT between the reactants ( 1 and 2a) in the TS configuration by around 30% (from −84.7 to −110.6 kJ/mol). To put it in perspective, this increase of −25.9 kJ/mol is more than three-fold larger than the aforementioned CT value in the isolated water trimer (−6.9 kJ/mol). The CT energy increase is associated with an increase of 61% of the total electron density transfer (from 184 to 297 millielectron, see Table S2 in the ESI).

The findings regarding the CT term are partially confirmed by the corresponding NBO calculations, see Figure S2 in the ESI. In agreement with ALMO EDA, the NBO analysis shows that the total CT energy between 2a and 3 is equivalent in all three configurations (i), (ii) and (iii). Nevertheless, according to the NBO results reported in Figure S2, there is only a marginal enhancement of CT between 1 and 2a, when the coordinates of water are included (an enhancement of 2.3 kcal/mol from −83.7 to −86.0 kcal/mol in the iv and ii configurations, respectively). This ambiguity can be attributed to the difference between CT in the NBO and ALMO EDA schemes. In NBO, the CT term is calculated from the off-diagonal blocks in the Fock-matrix of a given bond complex. This means that there are no variationally optimised intermediate states in the NBO scheme. Consequently, CT energy is not computed from a well-defined lowest-energy ‘zero-CT’ state as an accurate reference [9,21]. The ALMO constraint, however, strictly prohibits any CT between fragments at the optimised ${\mathrm{\Psi }}_{\mathrm{POL}}$ wavefunction, which is then considered to be the lowest-energy ‘zero-CT’ state [13]. Therefore, the quantity $\mathrm{\Delta }{E}_{\mathrm{CT}}$ in Equation 3(3) $\mathrm{\Delta }{E}_{\mathrm{CT}}=E({\mathrm{\Psi }}_{\mathrm{DEL}})-E({\mathrm{\Psi }}_{\mathrm{POL}})$(3) , accurately represents the energy of CT with no contribution from the polarisation energy [9,13,21].

Besides amplifying the CT within reactants, the 3 water molecules increase the total polarisation of the system from −29.6 to −70.8 kJ/mol. This increase should not be falsely attributed to stronger H-bonding in the TS. The polarisation term for the isolated interaction between 2a and 3 in iii is equivalent to that in i, see Table 1. In addition, there are favourable polarisation interactions within water molecules (−7.8 kJ/mol, see above) and between 1 and 3 (−5.0 kJ/mol, see Table 3). Hence, there seems to be no significant enhancement of the polarisation contribution to the extent of H-bonding in the TS as compared to the reactants.

5. Conclusions

Our results confirm that H-bonding to interfacial water molecules stabilises the TS of the cycloaddition of 1 and 2a. Nevertheless, factoring in all contributing terms to the strength of H-bonding, our EDA analysis shows that the H-bonds between the water molecules and 2a in both the reactant and TS configurations are of equivalent strengths, or even a bit stronger on the reactant side, contrary to what has been deduced by Jung and Marcus [4] The energy lowering originates from the enhanced CT between 1 and 2a and the more favourable polarisation and frozen density terms for the whole system of 1, 2a and 3. This work illustrates the use of EDA to investigate the role of interfacial water in on-water catalysis. Further work is needed to reach more decisive conclusions. The configurations used in EDA should be properly sampled from a simulation of heterogeneous water–organic system. Reactants that were predicted to yield inefficient on-water catalysis, like 2c in Figure 1 by Jung and Marcus [4], should be compared to 2b and 2a. Last but not least, one should also consider a reaction mechanism involving proton-transfer as a possible reason for the observed kinetic isotope effect in on-water catalysis.

Acknowledgements

This research has been enabled via funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 716142). In addition, we acknowledge the allocation of computing resources by the Paderborn Center for Parallel Computing (PC2).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research has been enabled via funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program [grant agreement number 716142].