Borane derivatives of five-membered N-heterocyclic rings as frustrated Lewis pairs: activation of CO2

The reaction of seventeen borane derivatives of five-membered N-heterocyclic rings (BNHRs) with CO2 has been studied by means of DFT calculations. Several non-covalent complexes between the BNHRs and CO2 which evolve through a TS in a single adduct for each BNHR have been identified. The calculated IRC of the TS has allowed to identify the non-covalent complex involved in the reaction in each case. The stationary points of the reactions have been analyzed with the distortion/interaction partition model. In addition, empirical models have been attempted to correlate the acid (fluoride ion affinity) and basic (proton affinity) properties of the isolated BNHR with the TS barriers and adduct energies. The energetics of the reactions are influenced by the number of nitrogen atoms in the ring.


Introduction
The continuous increase of carbon dioxide by human activities and its implication in the climate change [1,2] is generating challenging scientific opportunities to remove it from the atmosphere. Carbon capture and geological sequestration [3] is the only partial solution implemented nowadays at big scale. However, other alternatives have been proposed [4]. Among them is the capture and activation of carbon dioxide using different chemical reactions. Several types of compounds such as carbenes [5][6][7][8], phosphines [9][10][11], guanidines [12][13][14], and FLPs (frustrated Lewis pairs) [15][16][17] have been shown to interact with CO 2 and modify its geometrical/electronic properties making it suitable to react with other systems.
The reactivity and interaction of nitrogen heterocycles with CO 2 have been studied experimentally and theoretically. Thus, the formation of benzene-fused azole compounds through chemical fixation of CO 2 with aniline derivatives has been described [18], and the solubility of CO 2 in azole-based deep eutectic solvents has shown negative solution enthalpies [19,20], besides the possibility to use ionic liquids composed by imidazolium, pyridinium, and pyrrolidinium cations to fix CO 2 has been explored [21,22]. Regarding theoretical studies that have addressed the potential complexes between azoles and CO 2 , the complex between CO 2 and imidazole [23,24] shows three minima, two planar structures, and another with C 1 symmetry. A systematic study of the complexes of five-membered heterocycles with CO 2 [25][26][27][28][29] shows hydrogen, tetrel [30][31][32], and π-π interactions. The absorption of CO 2 in azole-based ionic liquids has been studied using a combination of DFT and molecular dynamics, with the formation of intermolecular N-C bonds between the nitrogen of the azole and the carbon of CO 2 [33]. Similar results are obtained when the interaction occurs between azole anions and CO 2 , leading to the activation of the latter [34,35].
In this article, seventeen borane derivatives of five-membered N-heterocyclic systems, with the borane group in α position to a pyridine-type nitrogen, have been studied as suitable molecules to activate CO 2 (Fig. 1). The nomenclature used to identify these compounds combines an identification of the 5-membered ring and the location of the BH 2 group. In some molecules, the interaction and reaction with CO 2 can occur in two different sides as indicated in Fig. 1 with the labels A and B. For each reaction, pre-reactive complexes and adducts energy minima have been characterized with the corresponding transition states (TS) connecting them through a reaction coordinate. A total of 20 reactions have been considered.

Computational methods
The geometry of the systems has been optimized with the M06-2X functional [49] and the aug'-cc-pVTZ basis set [50]. This basis set corresponds to the aug-cc-pVTZ for the heavy atoms and the cc-pVTZ one for the hydrogens. Frequency calculations at the same computational level have been carried out in the optimized geometries to confirm that they correspond to energy minima (zero imaginary frequencies) or transition states (only one imaginary frequency). The geometries of the stationary points have been gathered in Table S1 of the Supporting Information Material. The intrinsic reaction coordinate has been calculated to identify the minima structures connected through the located TS structures (Table S2). The distortion/interaction model has been used to partition the energy of the stationary points into the deformation of the monomers plus the interaction between the monomers in the geometry found at each stationary point [51]. All calculations were carried out with the Gaussian-16 program [52].
The molecular electrostatic potential (MESP) of the isolated azoles has been calculated and analyzed on the 0.001 au electron density isosurface with the Multiwfn program [53] and represented with the Jmol program [54]. Negative regions of MESP are associated with basic areas, while positive ones are associated with acid areas.
The acidity and basicity of the isolated FLP have been calculated by means of the proton affinity (PA) and fluoride ion affinity (FIA) [55], respectively. These magnitudes correspond to the negative values of the enthalpy change of the reactions 1 and 2, being in both cases positive values. In the case of the PA, the nitrogen involved in the reaction has been chosen for the protonation. The FIA always involves the BH 2 group.
The electron density of the systems has been analyzed within the quantum theory of atoms in molecules (QTAIM) [56,57] with the AIMAll program [58]. The points in space where the gradient of the electron density is null are defined as critical points (cp). Based on the sign of the curvatures, four types of cp can be defined as: nuclear attractors (3,-3), bond critical points, bcp, (3, -1), ring critical points (3, +1), and cage critical points (3, +3). This information complemented with the lines of the minimum gradient connecting the bcp with the nuclear attractors, bond paths, provides the molecular graph. In addition, the properties at the bcp give information on the characteristics of the corresponding bond or interaction.

Results and discussion
The initial exploration of the non-covalent complexes between the borane derivatives of five-membered N-heterocyclic rings (BNHR) with CO 2 provides up to three non-covalent complexes per BNHR molecule. In addition , two additional stationary points have been characterized in the reaction: the adduct and TS connecting one of the prereactive complexes with the adduct. The reaction of 1234-1H-5BH2 is an exception since an additional intermediate minimum and two TSs have been characterized.
The nomenclature used to differentiate the stationary points along the article will be BNHR:CO 2 for the noncovalent complexes, BNHR/CO 2 for the TSs, and BNHR-CO 2 for the adducts. The molecular graph of the three stationary points for the reaction of pz-1BH2 with CO 2 is shown in Fig. 2, and those for all the BNHR are included in Table S1 of the SI. This section has been divided into four parts. In the first part, the molecular and electronic characteristics of the isolated borane-azoles are analyzed. In the second part, the non-covalent complexes between the BNHR and CO 2 are discussed. The third part includes the analysis of the adducts and the last one the TSs of the reaction.

Isolated BNHR
The optimized structure of the isolated BNHR has C s symmetry with all the atoms in same plane. Ten of the studied BNHR show the borane group bonded to C atoms The borano-azoles considered in this work. Those with two potential reaction sides are indicated with A and B. Borano-azole labels are indicated underneath each structure with an average B-C distance of 1.53 Å, while in seven cases, a B-N bond is observed with an interatomic distance of 1.43 Å, an indication that the B-N bond is more polar and stronger than the B-C one [59]. In fact, comparing the energies of isomers of BNHR, those with B-N bonds are more stable than the corresponding B-C ones between 70 and 100 kJ mol −1 .
The BNHR presents a maximum of the MESP associated to the boron atom and a minimum on the nitrogen in α position (see two examples in Fig. 3 and the MESP of all the molecules in Fig. S1). The values of the MESP minima, σ min , range between −68 and −181 kJ mol −1 and those of the maxima, σ max , between +51 and +202 kJ mol −1 ( Table 1).
The larger values of the minima, in absolute value, are found in the molecules with a lower number of nitrogen atoms (imidazole and pyrazole) and the smallest in those molecules with more nitrogen atoms, like pentazole. The values of the maxima show the opposite trend; the largest is found in the pentazole derivative and the smallest in the imidazole. Acceptable linear correlations (R 2 > 0.91) are obtained between the minima and maxima values and the number of nitrogen atoms in the molecules with N-BH 2 bonds (Fig. S1).
The calculated fluoride ion affinity (FIA) and proton affinity (PA) of the isolated BNHR have also been gathered in Table 1. The FIA range between 261.7 and 420.3 kJ mol −1 and the PA between 932. 8

Non-covalent complexes
The potential non-covalent complexes between the BNHR and CO 2 in the region where the adduct formation can happen have been explored. Four different types of complexes have been found (Fig. 4): the first one corresponds to the tetrel bond interaction of the α-nitrogen to the borane group with the CO 2 atoms located approximately in the same plane. Depending on the groups in β position of the borane, additional hydrogen bonds with the CO 2 can be established in these complexes. In the second type of non-covalent complexes, the CO 2 is perpendicular to the plane of the BNHR with bifurcated tetrel bond between the Nα and one of the hydrogen atoms of the borane group. A similar situation is found in the third type of complexes; in those cases, there is a N in β to the borane group, and Nα and Nβ act simultaneously as tetrel bond acceptors [60]. The last case corresponds to a location of the CO 2 out of the BNMR plane interacting with the π-cloud or between one of the oxygen and the boron atom of the BNMR molecule. The most stable complex for each kind of interaction pattern is represented in Fig. 4.
The binding energies in the non-covalent complexes range between −9 and −22 kJ mol −1 ( Table 2). There is a modest correlation between the energy of the most stable noncovalent BNMR:CO 2 complexes and the number of nitrogen atoms in the BNMR; hence, this parameter is related to the binding energy, but the dispersion is still very large. The average values of the dissociation energies based on the BNHR with the same number of nitrogen atoms provide a good correlation with the number of N atoms (Fig. S5).

Adduct formation with CO 2
Once the non-covalent complex is formed, the CO 2 molecule can evolve and form an adduct with the BNHR through the corresponding TS. The relative energies of the adducts with Table 1 σ max associated to the BH 2 group, σ min of the N in α position to the BH 2 , fluoride ion affinity, FIA, and proton affinity, PA (kJ mol −1 ) calculated using Eqs. (3)    respect to the isolated BNHR and CO 2 have been gathered in Table 3. The adducts are more stable than the isolated BNHR and CO 2 in all cases except in 1234-2H-5BH2 (B) and 12345-1H-1BH2, the two latter presenting positive relative energies (+11.2 and +4.4 kJ mol −1 , respectively). The largest stabilization due to the formation of the adducts is found in the imi-2BH2 with a relative energy of −103 kJ mol −1 .
Interestingly, in the three most stable adducts (imi-2BH2, 124-4H-3BH2, and 124-1H-5BH2), the BH 2 is bonded to a carbon atom and presents a NH group in α' to the BH 2 , an amidine group surrounding the BH 2 . The remaining system with the same characteristics (1234-1H-5BH2) shows the 5th most stable adduct. The comparison of the relative energy of the adducts with the PA of the isolated BNHR provides a good linear relationship when the systems are divided in three groups: N-BH 2 molecules, C-BH 2 molecules with NH group in α', and the rest of the C-BH 2 molecules (Fig. S6). No significant linear relationship has been found between the relative energy and the FIA of the isolated BNHR. However, when the PA and FIA were considered together (Eqs. 3 and 4), an acceptable multiple relationship with R 2 = 0.96 is obtained for all the systems simultaneously (Fig. 5). In order to elucidate the importance of FIA and PA, their The distortion/interaction partition terms (Table 3) show very large interaction energies between the two molecules (between −319 and −499 kJ mol −1 ) that are compensated by the distortion of the CO 2 (between +190 and +232 kJ mol −1 ) and the azole (between +138 and +201 kJ mol −1 ). The representation of the distortion energies vs. the relative energy of the adducts in the systems with C-BH 2 and N-BH 2 bonds (Fig. 6) shows that  in each of these two families, the distortion of the CO 2 is larger than that of the azole. In addition, the distortion values for the systems with C-BH 2 bonds are smaller than the ones with N-BH 2 with similar relative energy of the adducts. Good linear correlations are found between these parameters (R 2 > 0.96) except for the distortion of the azoles in the systems with C-BH 2 bond and the relative energy of the adduct (R 2 = 0.38).

TS of the adduct formation
The energy profile of the CO 2 addition to the BNHR shows a unique TS connecting one of the non-covalent complexes with the adduct except for the reaction of 1234-1H-5BH2 that presents two TS and an intermediate linking them (Fig. 7). The intermediate in this reaction is only 0.6 kJ mol −1 more stable than TS2. The relative energies of the TS that connect the noncovalent complexes with the adducts range between −15 and +31 kJ mol −1 (Table 4). In three cases, the relative energy of the TS is negative (124-1H-5BH2, 124-4H-3BH2, and 1234-1H-5BH2) which indicates that they are more stable than the isolated BNHR and CO 2 . The calculation of the IRC profiles allows to identify the non-covalent complex directly involved in the adduct formation through the TS (in bold in Table 2). In eleven cases, they correspond to the "out of plane" complexes, in seven to the Nα complex, and in two cases to the Nα-BH one. A previous study [48] that does not consider the IRC calculations assumes that the prereactive complex for the imidazole and pyrazole derivatives is the Nα complex, while our calculations show alternative pre-reactive complexes for the pyrazole derivatives. The intermolecular N···C distances in the TS range from 1.90 and 2.66 Å and the B···O between 1.61 and 2.29 Å in the TS structures. Reasonable linear relationships (Fig. S7) are found between these geometrical parameters with positive slope, an indication that both distances increase simultaneously.
The analysis of the N···C and B···O interatomic distances in the TS clearly indicates that the reaction proceeds asynchronously since in the TS the B-O bonds are shorter than the N-C ones, between 0.69 and 0.24 Å. Comparing the interatomic distances in the TS with those of the adducts, it is observed that the B-O distances in the TS are on average 1.13 times larger than the corresponding ones in the adduct, while the analogous relationship in the C-N bond is 1.52.
Thus, the formation of the B-O bond is more advanced in the TS than the C-N one.
The AIM analysis of the electron density at the TS shows intermolecular bcp for the N···C and B···O contacts. The values of the electron density properties at the N···C bcp are typical of interactions that start to show incipient covalent character with ρ bcp between 0.09 and 0.02 au, positive values of the Laplacian and small values of H, being positive for the longest distances and negative for the shortest ones.
The B···O contacts show ρ bcp values between 0.09 and 0.04 au, with positive values of the Laplacian due to the large electronegativity of the atoms involved (it is also positive in the short B-O bonds found in the adducts) but negative values of H for all the cases, which confirms that this bond is more formed in the TS than the N···C one.
The distortion/interaction partition of the systems in the TS (Table 4) shows that the distortion of the azole is larger than the one of the CO 2 in contrast to what is observed in the adducts.

Conclusions
The reaction of seventeen borane derivatives of fivemembered N-heterocyclic rings (BNHRs) with CO 2 , which yield 20 adducts, has been studied by means of DFT methods. The acid/base properties of the isolated BNHR compounds have been characterized by means of the fluoride ion affinity (FIA) and proton affinity (PA). These properties are related to the values obtained for the molecular electrostatic potential maxima and minima in the corresponding atomic centers, respectively.
The reaction evolves in 19 of the 20 cases studied with the formation of a non-covalent complex, a TS, and an adduct. The remaining case (1234-1H-5BH2) presents an additional intermediate and two TSs connecting the non-covalent complex with the adduct. Four different types of non-covalent complexes have been characterized between the BNHRs and CO 2 with binding energies up to −20 kJ mol −1 . The IRC calculations have clarified the non-covalent complex directly involved in the reaction. The TS show relative energies up to +31 kJ mol −1 with several cases with negative values. The adducts are more stable than the isolated BNHR and CO 2 except in three cases. The most favorable cases present an NH group in α' of the C-BH 2 moiety with relative energies with respect to the entrance channel up to −103 kJ mol −1 .
Overall, the stability of the adducts and the small values of the TSs indicate that BNHRs, based on the systems studied here, could be used to sequestrate and activate CO 2 molecules, especially in cases with two or three nitrogen atoms since the increment of nitrogen atoms disfavors the stability of the adducts.

Acknowledgements
This work was carried out with financial support from the Ministerio de Ciencia, Innovación y Universidades (Project PID2021-125207NB-C32) and Comunidad de Madrid (PS2018/EMT-4329 AIRTEC-CM). Thanks are given to the CTI (CSIC) and to the Irish Centre for High-End Computing (ICHEC, Dublin) for their continued computational support.
Author contribution M. Ferrer and I. Alkorta performed the calculations. I. Alkorta wrote the first draft. All the authors revised the manuscript.
Availability of data and materials All data generated or analyzed during this study are included in this published article and its supplementary information files.

Declarations
Ethical approval Not applicable.