Non-Covalent Interactions in Hydrogen Storage Materials LiN(CH₃)₂BH₃ and KN(CH₃)₂BH₃

Abstract

In the present work, an in-depth, qualitative and quantitative description of non-covalent interactions in the hydrogen storage materials LiN(CH₃)₂BH₃ and KN(CH₃)₂BH₃ was performed by means of the charge and energy decomposition method (ETS-NOCV) as well as the Interacting Quantum Atoms (IQA) approach. It was determined that both crystals are stabilized by electrostatically dominated intra- and intermolecular M∙∙∙H–B interactions (M = Li, K). For LiN(CH₃)₂BH₃ the intramolecular charge transfer appeared (B–H→Li) to be more pronounced compared with the corresponding intermolecular contribution. We clarified for the first time, based on the ETS-NOCV and IQA methods, that homopolar BH∙∙∙HB interactions in LiN(CH₃)₂BH₃ can be considered as destabilizing (due to the dominance of repulsion caused by negatively charged borane units), despite the fact that some charge delocalization within BH∙∙∙HB contacts is enforced (which explains H∙∙∙H bond critical points found from the QTAIM method). Interestingly, quite similar (to BH∙∙∙HB) intermolecular homopolar dihydrogen bonds CH∙∙∙HC appared to significantly stabilize both crystals—the ETS-NOCV scheme allowed us to conclude that CH∙∙∙HC interactions are dispersion dominated, however, the electrostatic and σ/σ*(C–H) charge transfer contributions are also important. These interactions appeared to be more pronounced in KN(CH₃)₂BH₃ compared with LiN(CH₃)₂BH₃.

1. Introduction

An increase in energy consumption as well as the environmental harmfulness of current coal or hydrocarbon based fuels has led to intensive search for alternative energy sources [1,2,3]. Therefore, various hydrogen storage materials, that contain significant amounts of hydrogen, have been recently proposed. Boranes are probably one of the best known group among numerous hydrogen storage materials [4,5,6,7,8,9,10]. For example one can present ammonia borane (NH₃BH₃) [4,5,6,7,8,9,10]—the attractiveness of this material stems from its high stability, even at higher temperature (the melting point is 104 °C), as well as its large hydrogen storage capacity (19.6 wt% H₂). It has been demonstrated that the former feature of ammonia borane crystal originates predominantly from the existence of polar dihydrogen bonds N–H^δ+∙∙∙^−δH–B between monomers [11,12,13,14,15,16]. Furthermore, it has been proven that the presence of N–H^δ+∙∙∙^−δH–B as well as other non-covalent interactions not only determines the stability, but it can also facilitate various steps of dehydrogenation [5,6,7,8,11,12,13,14,15,16,17,18,19,20,21,22,23].

It has been shown that incorporation of alkali metals into boranes might accelerate thermolitic dehydrogenation as well as reduce formation of volatile byproducts [24,25,26]. Therefore various hybrid type materials have also been proposed and extensively studied—as examples one can present LiBH₄/NH₃BH₃ [27], M[Zn(BH₄)₃], M = Li, Na, K [28] or Al(BH₄)₃·NH₃BH₃ [25].

Quite recently McGrady and coworkers published the cutting-edge article in which they had synthesized and characterized two further hydrogen rich crystals LiN(CH₃)₂BH₃ and KN(CH₃)₂BH₃ [29]. They are depicted in Figure 1. In addition, the authors performed topological electron density based study by means of the QTAIM method of Bader [30]—it was reported that mainly M∙∙∙H–B (M = Li, K) interactions stabilize the crystals. Remarkably, the authors also emphasized that, apart from the mentioned non-covalent interactions, one observes untypical homopolar dihydrogen interactions of the types BH∙∙∙HB and CH∙∙∙HC which are found to determine the chain-like 1D architecture of LiN(CH₃)₂BH₃ and 2D layers in KN(CH₃)₂BH₃ crystal [29], Figure 1 and Figure S1. It is noteworthy that these types of connections are intuitively considered as destabilizing due to the lack of electrostatic attraction between hydrogen atoms—homopolar dihydrogen interactions (especially the intramolecular ones) are still a matter of debate in the literature [31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46]. Very recently more and more evidence has been reported in the literature that highlight the stabilizing nature of homopolar dihydrogen interactions [17,21,29,36,37,38,39,40,41,42,43,44,45,46].

Accordingly, in this work we provide complementary results which shed light on energetic, quantitative and qualitative characteristics of non-covalent interactions that contribute to the stability of LiN(CH₃)₂BH₃ and KN(CH₃)₂BH₃ crystals. It is an important goal as it is known that purely topological QTAIM analysis might provide bond paths between atoms (or fragments) even in situations where the overall interaction energies are positive (destabilizing) [36]. Therefore, we applied the charge and energy decomposition scheme (ETS-NOCV) [47,48,49] which has been proven to provide compact, qualitative and quantitative, descriptions of various types of chemical bonds starting from strong covalent bonds, going through dative connections [47,48,49] and ending up at various non-covalent interactions [7,8,37,50,51]. We applied the ADF program [52,53,54] in which the ETS-NOCV scheme is implemented. In order to achieve our goal we chose two types of cluster models—the first type (containing four monomers), marked by a blue dashed line in Figure 1A, is suitable for extraction of M∙∙∙H–B (M = Li, K), and BH∙∙∙HB interactions, whereas the second one (containing eight monomers), depicted by a red dashed line (Figure 1B), contains the CH∙∙∙HC contacts. For selected models we also plotted the reduced density gradient of the NCI (Non-Covalent Index) method [55] in order to qualitatively characterize non-covalent interactions. In addition, the Interacting Quantum Atoms (IQA) energy decomposition scheme [56] was applied for a quantitative description of selected non-covalent interactions in LiNMe₂BH₃.

2. Materials and Methods

Our calculations were performed by means of the Amsterdam Density Functional (ADF) program [52,53,54]. We used DFT/BLYP-D3/TZP because it has been proven many times in the past that they provide satisfactory results for non-covalent interactions [57,58,59]. The empirical Grimme correction (D3) [60] was used as implemented in the ADF program. We did not calculate basis set superposition errors (BSSE) because these effects are captured in the empirical correction D3 [60]. It is important to emphasize that we also performed additional test calculations for the tetrameric cluster of LiN(CH₃)₂BH₃ using the following methods: PBE-D3/TZP, BP86-D3/TZP (ADF program) as well as MP2/6-311 + G**, PBE-D3/6-311 + G**, BP86-D3/6-311 + G**, MO6-2X/6-311 + G**, wB97XD/6-311 + G** (based on the Gaussian package) [61]—the results of total bonding energies appeared to be very similar among all these methods, Tables S1,S2 in ESI file. Accordingly, in the main text we have only discussed the results from ADF/DFT/BLYP-D3/TZP.

In the next paragraph the main formulas of charge and energy decomposition scheme ETS-NOCV are outlined.

2.1. ETS-NOCV

ETS-NOCV method [49] is a merger of the energy decomposition scheme by Ziegler-Rauk [53,62] with the Natural Orbitals for Chemical Valence (NOCV) [47,48].

The natural orbitals for chemical valence (NOCV) are eigenvectors that diagonalize the deformation density matrix:

$$\mathsf{\Delta }P{C}_i=v_i{C}_i,{\mathsf{\Psi }}_i={\displaystyle {\sum }_j^N{C}_{i,j}{\lambda }_j}$$

(1)

where C_i is a vector of coefficients, expanding Ψ_i in the basis of fragment orbitals λ_j; N is a total number of fragment λ_j orbitals. It was shown that the natural orbitals for chemical valence pairs (ψ_–k,ψ_k) decompose the differential density Δρ into NOCV-contributions (Δρ_k):

$$\mathsf{\Delta }\rho (r)={\displaystyle \sum _{k=1}^{M/2}{v}_k\left[-{\psi }_{-k}^2(r)+{\psi }_k^2(r)\right]={\displaystyle \sum _{k=1}^{M/2}\mathsf{\Delta }{\rho }_k(r)}}$$

(2)

where v_k and M stand for the NOCV eigenvalues and the number of basis functions, respectively. Visual inspection of deformation density plots (Δρ_k) helps to attribute symmetry and the direction of the charge flow. In addition, these pictures are enriched by providing the energetic estimations, ΔE_orb(k), for each Δρ_k within the ETS-NOCV scheme.

The exact formula, which links the ETS and NOCV methods, is given in the next paragraph, after we briefly present the basic concept of the ETS scheme. In this method the total bonding energy, ΔE_total, between interacting fragments is divided into four components:

ΔE_total = ΔE_dist + ΔE_elstat + ΔE_Pauli + ΔE_orb = ΔE_dist + ΔE_int

(3)

One could add that the negative value –ΔE_total is a bond dissociation energy. The first contribution of Equation (3), ΔE_dist, describes the amount of energy required to promote fragments from their equilibrium geometry to the conformations they adopt in the final optimized molecule. The second term, ΔE_elstat, corresponds to the classical electrostatic interaction between the promoted fragments as they are brought to their positions in the final complex. The third term, ΔE_Pauli, accounts for the repulsive Pauli interaction between occupied orbitals on the two fragments in the combined molecule. Finally, the last stabilizing term, ΔE_orb, represents interactions between the occupied molecular orbitals of one fragment with the unoccupied molecular orbitals of the other fragment as well as mixing of occupied and virtual orbitals within the same fragment (inner-fragment polarization). This energy term is linked to the electronic bonding effect coming from the formation of a chemical bond. In the combined ETS-NOCV scheme [49] the orbital interaction term (ΔE_orb) is expressed in terms of NOCV’s eigenvalues (v_k) as:

$$\mathsf{\Delta }{E}_{orb}={\displaystyle \sum _k\mathsf{\Delta }{E}_{orb}(k)={\displaystyle \sum _{k=1}^{M/2}{v}_k\left[-F_{-k,-k}^{TS}+F_{k,k}^{TS}\right]}}$$

(4)

where $F_{i,i}^{TS}$ are diagonal Kohn-Sham matrix elements defined over NOCV with respect to the transition state (TS) density at the midpoint between the density of the molecule and the sum of fragment densities. The above components ΔE_orb(k) provide the energetic estimation of Δρ_k that may be related to the importance of a particular electron flow channel for the bonding between the considered molecular fragments. Finally, in this work we applied a dispersion corrected functional, so this term (ΔE_disp) enters additionally into Equation (3). The ETS-NOCV analysis was done based on the Amsterdam Density Functional (ADF) package in which this scheme was implemented.

For ETS-NOCV analyses the crystal coordinates (not reoptimized) have been predominantly used in order to reflect the real structures of LiN(CH₃)₂BH₃ and KN(CH₃)₂BH₃. Accordingly, the distortion energy term, ΔE_dist (Equation (3)), was not calculated. Hence, an interaction energy ΔE_int was mostly analyzed in this work. We found for the tetrameric lithium model that the reoptimized geometry as well as bonding properties corresponds well to the crystal structure—accordingly, for this system we also considered the ΔE_dist term (Figure 2A) and discussed the corresponding total bonding ΔE_total and interaction ΔE_int energies. It needs to be added that our efforts to reoptimize the remaining models provided geometries which do not correspond to the crystal structures of LiN(CH₃)₂BH₃ and KN(CH₃)₂BH₃. For monomers the optimized structures are considered for ETS-NOCV analysis as they are similar to crystal geometries—in addition, it allows for discussion of the overall stability of monomers (ΔE_total values).

2.2. NCI Technique

It has been shown that the reduced density gradient:

$$s=\frac{1}{2{(3{\pi }^2)}^{1/3}}\frac{\left|\nabla \rho \right|}{{\rho }^{4/3}}$$

(5)

appeared to be a useful quantity for a description of non-covalent interactions. In order to obtain information about the type of bonding, plot of reduced density gradient s against molecular density ρ is very often examined. When a weak inter- or intramolecular interaction is present, there exists a characteristic spike lying at low values of both density ρ and reduced density gradient s. To distinguish between attractive and repulsive interactions the eigenvalues (λi) of the second derivative of density (Hessian, ∇²ρ) are used, ∇²ρ = λ1 + λ2 + λ3. Namely, bonding interactions are characterized by λ2 < 0, whereas λ2 > 0 indicates that the atoms are in non-bonded contact. Therefore, within the NCI technique, one can draw information about non-covalent interactions from the plots of sign(λ2)ρ vs. s. In such plots the low gradient spike, an indicator of stabilizing interaction, is located within the region of negative values of the density. On the contrary, the repulsive interactions are characterized by positive values of sign(λ2)ρ. One can also plot the contour of s colored by the value of sign(λ2)ρ providing a pictorial representation of non-covalent interactions.

2.3. IQA (Interacting Quantum Atoms) Energy Decomposition Scheme

The Interacting Quantum Atoms (IQA) approach of Blanco and coworkers [56] allows to partition an electronic energy E into atomic (${\mathrm{E}}_{\text{self}}^{\mathrm{A}}$) and diatomic contributions (${\mathrm{E}}_{\text{int}}^{\text{AB}}$):

$$\mathrm{E}={\displaystyle \sum _{\mathrm{A}}{\mathrm{E}}_{\text{self}}^{\mathrm{A}}}+\frac{1}{2}{\displaystyle \sum _{\mathrm{A}}{\displaystyle \sum _{\mathrm{B}\ne \mathrm{A}}{\mathrm{E}}_{\text{int}}^{\text{AB}}}}$$

(6)

The interatomic interaction energy ${\mathrm{E}}_{\text{int}}^{\text{AB}}$ covers all inter-particle interactions: nucleus-nucleus, ${\mathrm{V}}_{\text{nn}}^{\text{AB}}$, nucleus-electron, ${\mathrm{V}}_{\text{ne}}^{\text{AB}}$, electron-nucleus, ${\mathrm{V}}_{\text{en}}^{\text{AB}}$, and electron-electron, ${\mathrm{V}}_{\text{ee}}^{\text{AB}}$, coming from interatomic interaction energies of particles ascribed to atom A with particles ascribed to atom B:

$${\mathrm{E}}_{\text{int}}^{\text{AB}}={\mathrm{V}}_{\text{nn}}^{\text{AB}}+{\mathrm{V}}_{\text{en}}^{\text{AB}}+{\mathrm{V}}_{\text{ne}}^{\text{AB}}+{\mathrm{V}}_{\text{ee}}^{\text{AB}}={\mathrm{V}}_{\text{nn}}^{\text{AB}}+{\mathrm{V}}_{\text{en}}^{\text{AB}}+{\mathrm{V}}_{\text{ne}}^{\text{AB}}+{\mathrm{V}}_{\text{eeC}}^{\text{AB}}+{\mathrm{V}}_{\text{eeX}}^{\text{AB}}$$

(7)

The ${\mathrm{V}}_{\text{ee}}^{\text{AB}}$ term can be further divided into exchange (${\mathrm{V}}_{\text{eeX}}^{\text{AB}}$) and Coulomb (${\mathrm{V}}_{\text{eeC}}^{\text{AB}}$) contributions. The AIMALL program was used for the IQA calculations [63]. Due to the fact that we are interested in interaction energies in crystals we focused our attention on ${\mathrm{E}}_{\text{int}}^{\text{AB}}$. More details and numerous applications of the IQA technique can be found in Reference [56].

2.4. Molecular Electrostatic Potential (MEP)

The electrostatic potential V(r) of a molecule at point “r”, due to nuclei and electrons, is defined as:

$$V(r)={\displaystyle \sum _A\frac{Z_A}{\left|R_A-r\right|}-{\displaystyle \int \frac{\rho (r\prime )dr\prime }{\left|r-r\prime \right|}}}$$

(8)

where Z_A is the charge of nucleus at position R_A and ρ (r) is the total electronic density. The sign of V(r) depends upon whether the positive contribution of the nuclei or negative from the electrons is dominant. Negative values of V(r) correspond to nucleophilic areas of the molecule, whereas the positive to electrophilic regions. It has been demonstrated in numerous works that MEP is a very useful quantity for in depth description of electron density distribution [64,65,66,67].

3. Results and Discussion

Let us start by discussing factors determining the stability of LiN(CH₃)₂BH₃. It can be seen from Figure 2A that the overall interaction energy between fragments, each consisting of the two LiN(CH₃)₂BH₃ monomers, is ΔE_int = −35.63 kcal/mol (−17.81 kcal/mol permonomer-monomer interaction). This is in very good agreement with ΔE_int obtained for the crystal (non-optimized) geometry, Figure S2. An inclusion of the geometry distortion term leads to an overall bonding energy, ΔE_total = −29.57 kcal/mol. It is noteworthy that other computational protocols provide very similar ΔE_total values, Table S1 and Table S2.

The dominating contribution (55%) to ΔE_total stems from the electrostatic stabilization ΔE_elstat = −35.56 kcal/mol, followed by the orbital interaction ΔE_orb = −18.17 kcal/mol (28%) and the dispersion components ΔE_disp = −11.38 kcal/mol (17%), Figure 2A. The molecular electrostatic potential (MEP) of the monomer demonstrates (Figure 3A) that the borane group is negatively charged, whereas the electrophilic region (positive MEP values) is seen around Li which explains the dominance of the electrostatic term in the intermolecular stabilization of LiN(CH₃)₂BH₃.

Further decomposition of ΔE_orb into deformation density contributions according to the ETS-NOCV scheme is shown in Figure 2B. It can be seen that the leading deformation densities, Δρ₁ and Δ₂, depict the formation of Li∙∙∙H–B bonds. They originate from an outflow of electron density from σ (B–H) bonds to Li⁺ ions, Figure 2B, and correspond to the following stabilization, ΔE_orb(1) = −11.91 kcal/mol, ΔE_orb(2) = −3.47 kcal/mol, respectively. Such outflow leads to the elongation of B-H bonds by ~0.1 Å (as compared with non-bonding monomers). It is important to emphasize that Li∙∙∙H–B interactions also lead to some charge delocalization within the “bay” containing two Li ions and two borane units, see Δρ₂ in Figure 2B. This is fully consistent with the presence of QTAIM bond critical points between hydrogen atoms in homopolar bridges BH∙∙∙HB as found by McGrady and coworkers [29]. However, an interesting question emerges at this point: is the LiN(CH₃)₂BH₃ dimer consisting of two monomers exposed to each other via pure BH∙∙∙HB interaction stable? Such a situation, i.e., the dimer in the geometry of the optimized tetrameric cluster model, is depicted in Figure 4B. Our results clearly indicate that in such a case the overall monomer-monomer interaction energy is positive due to significant Pauli and electrostatic repulsions (with total of +9.06 kcal/mol) that overcome subtle stabilization from charge transfer ΔE_orb and dispersion ΔE_disp (summing up to −4.16 kcal/mol)—accordingly, pure BH∙∙∙HB interactions would destabilize the LiN(CH₃)₂BH₃ crystal. The same conclusions are valid when considering the dimer in crystal geometry, Figure S3.

Furthermore, we have not found any local minimum energy structure on the potential energy surface for a dimer that would contain solely BH∙∙∙HB interaction (a dimer where BH₃ units are exposed to each other). Summarizing, the major stabilization in the LiN(CH₃)₂BH₃ crystal arises from electrostatically dominated intermolecular Li∙∙∙H–B interactions, as previously reported by McGrady et al. [29]. The formation of such bonds additionally enforce some charge delocalization within BH∙∙∙HB contacts (Δρ₂ in Figure 2B), however, the overall BH∙∙∙HB interaction is destabilizing. Our conclusions are fully in line with the other, mostly experimental studies, in which the destabilizing role of BH∙∙∙HB interactions were also suggested [68,69,70,71,72]. McGrady and coworkers reported the opposite in their series of recent articles [17,20,21,29]. Due to the fact that such interactions are clearly a matter of debate, in addition, we performed the Interacting Quatum Atoms (IQA) based study for the tetrameric models of LiNMe₂BH₃. The results are gathered in

Table 1. The interacting quatum atoms (IQA) energy decomposition results (in kcal/mol) describing the two atomic interactions X∙∙∙Y (X=Li, H; Y=H) in LiNMe₂BH₃.

IQA(X∙∙∙Y)	${\mathrm{V}}_{\text{ne}}^{\text{AB}}$	${\mathrm{V}}_{\text{en}}^{\text{AB}}$	${\mathrm{V}}_{\text{nn}}^{\text{AB}}$	${\mathrm{V}}_{\text{ee}}^{\text{AB}}$	${\mathrm{V}}_{\text{eeC}}^{\text{AB}}$	${\mathrm{V}}_{\text{eeX}}^{\text{AB}}$	${\mathrm{E}}_{\text{int}}^{\text{AB}}$*
Li∙∙∙H(B)	–805.9	–338.4	480.1	565.4	568.2	–2.8	–98.8
(C)H∙∙∙H(C)	–140.6	–138.2	130.8	147.2	148.6	–1.4	–0.80
(B)H∙∙∙H(B)	–198.9	–198.6	119.7	327.2	330.0	–2.7	+49.4

* ${\mathrm{E}}_{\text{int}}^{\text{AB}}={\mathrm{V}}_{\text{nn}}^{\text{AB}}+{\mathrm{V}}_{\text{en}}^{\text{AB}}+{\mathrm{V}}_{\text{ne}}^{\text{AB}}+{\mathrm{V}}_{\text{ee}}^{\text{AB}}={\mathrm{V}}_{\text{nn}}^{\text{AB}}+{\mathrm{V}}_{\text{en}}^{\text{AB}}+{\mathrm{V}}_{\text{ne}}^{\text{AB}}+{\mathrm{V}}_{\text{eeC}}^{\text{AB}}+{\mathrm{V}}_{\text{eeX}}^{\text{AB}}$.

.

It is visible from Table 1 that the overall diatomic interaction energies, ${\mathrm{E}}_{\text{int}}^{\text{AB}}$, for Li∙∙∙H(B) and (C)H∙∙∙H(C) are negative, ${\mathrm{E}}_{\text{int}}^{\text{Li}\cdots \text{H(B)}}=-98.8\text{ kcal}/\text{mol}$, ${\mathrm{E}}_{\text{int}}^{\text{(C)H}\cdots \text{H(C)}}=-0.80\text{ kcal}/\text{mol}$, which indicates the stabilizing interactions as opposed to the homopolar (B)H∙∙∙H(B) contacts which appeared to be significantly destabilizing, ${\mathrm{E}}_{\text{int}}^{\text{(B)H}\cdots \text{H(B)}}=+49.4\text{ kcal}/\text{mol}$. This is due to the significantly positive electron-electron repulsion term, ${\mathrm{V}}_{\text{ee}}^{\text{AB}}=327.2\text{ kcal}/\text{mol}$, caused in turn by the Coulomb contribution ${\mathrm{V}}_{\text{eeC}}^{\text{AB}}=330.0\text{ kcal}/\text{mol}$, Table 1. We also looked at the partial charges in the monomer and tetramer of LiN(CH₃)₂BH₃ and found that in both cases the borane units are negatively charged, which conforms to the significant value of the Coulomb repulsion ${\mathrm{V}}_{\text{eeC}}^{\text{AB}}$ found from the IQA analysis, Figure S4 and Table 1. We further calculated the electrostatic interaction between the monomers bonded via BH∙∙∙HB contacts in the presence of the two other monomers (Figure S5). Interestingly, we found the electrostatic repulsion, ca. 15.6 kcal/mol, which is even more pronounced than the repulsion noted for the dimer without the neighboring monomers, 7.08 kcal/mol, Figures S3,S5—we determined that it is due to the closely located Li ions which act as electron density attractors making BH₃ units more nucleophilic compared with the monomer, Figure S6. All these results based on the ETS-NOCV, IQA methods, atomic charges and molecular electrostatic potentials allow to conclude on the destabilizing nature of BH∙∙∙HB contacts in LiN(CH₃)₂BH₃—it is important to admit that our results are based on cluster models which might lead to omission of some bonding features in real crystals especially as far as weak non-covalent interactions are taken into account—accordingly, further studies based on other methods within the periodic calculations are necessary in order to fully delineate the nature of homopolar BH∙∙∙HB interactions. Unfortunately, the ETS-NOCV and IQA schemes are not yet available for public use in periodic calculations codes. On the other hand, it seems rational to comment, based on the very huge positive value of the (B)H∙∙∙H(B) interaction energy ${\mathrm{E}}_{\text{int}}^{\text{(B)H}\cdots \text{H(B)}}=+49.4\text{ kcal}/\text{mol}$ in the tetrameric Li-model, that it is to be expected that in the real Li-crystal the BH∙∙∙HB interactions are likely to be destabilizing.

It is probably important to add that the electron density accumulation in the inter-hydrogen region of BH∙∙∙HB is indeed sufficient to observe a bond path (as it has been noted by McGrady et al. [29])—however, it does not necessarily imply the overall stabilizing interactions: for example, as nicely demonstrated by Cukrowski et al. [36], two water molecules enforced to approach each other via oxygen atoms leads to formation of the oxygen-oxygen bond critical point, however, the overall binding energy is as expected positive (destabilizing) due to the fact that some subtle stabilization arising from the electron-exchange channel (resulting also in the negative ΔE_orb values) is diminished by the Pauli and electrostatic repulsion. A similar situation is observed in dimers of M₂X₂ (for M = Li, K, X = H, Cl) where M-X stabilization outweighs the repulsion stemming from M-M and X-X interactions [73].

In Figure 3A the structure of the LiN(CH₃)₂BH₃ monomer is presented. It can be seen that lithium ion forms a chemical bond not only with the nitrogen atom, but also with the BH₃ unit through intramolecular Li∙∙∙H–B interactions. The binding energy of BH₃ to NMe₂Li appears to be significant, ΔE_total = −63.99 kcal/mol, Figure 3B. Just for comparison, the binding energy of BH₃ to ammonia in NH₃BH₃ is only −30.3 kcal/mol [13]. Such a difference is related to the existence of strong intramolecular Li∙∙∙H–B interactions. ETS-NOCV allowed us to conclude that, apart from strong dative bonds, described by Δρ(N-B) and the corresponding ΔE_orb(B–N) = −105.63 kcal/mol, additional intramolecular Li∙∙∙H–B interactions are formed, Δρ(Li∙∙∙H–B). The latter component corresponds to significant charge transfer stabilization, ΔE_orb(Li∙∙∙H–B) = −20.94 kcal/mol, Figure 3B. It is noteworthy that the intramolecular Li∙∙∙H–B charge transfer component (Figure 3B) appears to be stronger than the intermolecular one, ΔE_orb(Li∙∙∙H–B) = −15.38 kcal/mol (Figure 2B)—this is related to the fact that in the latter case B–H bonds interact with the single Li ion, whereas in the former one with multiple Li ions.

We further performed similar calculations for the tetrameric model of KN(CH₃)₂BH₃ as well as for the monomer—the results of ETS-NOCV calculations are shown in Figure 5 and Figure 6 . Similarly to LiN(CH₃)₂BH₃, intermolecular interactions in the tetramer of KN(CH₃)₂BH₃ are dominated by electrostatic forces (60% of total stabilizing contributions), followed by dispersion (22%) and orbital interaction components (18%), Figure 5. The overall intermolecular K∙∙∙H–B interaction energy is slightly less pronounced compared to the corresponding Li∙∙∙H–B due to the larger size of potassium compared to lithium. One should point out the lack of BH∙∙∙HB electron delocalization within the “bay” between the two potassium atoms and BH₃ groups, which is present in LiN(CH₃)₂BH₃ (Δρ₂ in Figure 2b). It has been already noted by McGrady and coworkers [29]. The B–N dative bond in KN(CH₃)₂BH₃ monomer is of similar strength ΔE_orb(B–N) = −104.30 kcal/mol with respect to LiN(CH₃)₂BH₃. The intramolecular K∙∙∙H–B charge delocalization is of similar magnitude to the intermolecular-one, Figure 5B and Figure 6B. Finally, the overall bonding energy of BH₃ to KN(CH₃)₂ in the monomer, ΔE_total = −67.34 kcal/mol, appeared to be more negative compared to the corresponding ΔE_total = −63.99 kcal/mol for LiN(CH₃)₂BH₃ predominantly due to a lower Pauli repulsion contribution, Figure 3 and Figure 6.

Finally, we consider the two remaining cluster models of LiN(CH₃)₂BH₃ and KN(CH₃)₂BH₃ where homopolar CH∙∙∙HC interactions are involved, Figure 1B, Figure 7A and Figure 8A.

The results of energy decomposition analyses demonstrate the stabilizing nature of CH∙∙∙HC interactions in both crystals—the overall interaction energy is ΔE_int = −17.45 kcal/mol for KN(CH₃)₂BH₃, whereas it is only ΔE_int = −4.34 kcal/mol for LiN(CH₃)₂BH₃, Figure 7A and Figure 8A. Significantly more pronounced stabilization in potassium crystal originates from larger number of CH∙∙∙HC contacts as compared to the lithium analog. It constitutes the 2D layers architecture of the potassium crystal as compared to rather 1D chain-like structure in LiN(CH₃)₂BH₃. In both cases the main contribution (64%–68%) to overall stabilization is dispersion, ΔE_disp = −4.63 kcal/mol for LiN(CH₃)₂BH₃ and ΔE_disp = −18.78 kcal/mol for KN(CH₃)₂BH₃, Figure 7A and Figure 8A. The same has been already suggested by McGrady and coworkers [29]. We found in addition that the electrostatic ΔE_elstat and orbital interaction ΔE_orb terms are also important, Figure 7A and Figure 8A. Figure 7B and Figure 8B shows that the formation of the CH∙∙∙HC interactions is accompanied by a charge outflow from the occupied σ(C–H) bonds and the electron density accumulation is visible in the inter-hydrogen region of CH∙∙∙HC units. Considering the classical language of a molecular orbital theory one can summarize that the “electronic” part of the CH∙∙∙HC bonding is based on both donation from the occupied σ (C–H) bonds into the empty σ*(C–H) of methyl groups as well as polarization of the C–H bonds (mixing of σ/σ*(C–H)). These stabilizing contributions (the polarization and charge transfer) are clearly mixed—at this point one can reference other important and interesting works in which the meaning of both contributions is debated in terms of non-covalent interactions [67,68,74,75,76].

It is important to reference other works where dispersion dominated CH∙∙∙HC interactions have been found—Echeverría and coworkers found, based on MP2 studies, the “subtle but not faint” stabilizing CH∙∙∙HC interactions between alkanes [77]. Further Valence Bond studies have revealed [78], in line with our conclusions, that, apart from the crucial dispersion term, also the σ/σ*(C–H) polarization/charge transfer and electrostatic contributions are important. Recently, numerous reports are present in the literature on the stabilizing nature of CH∙∙∙HC interactions in various types of hydrocarbons [37,38,79,80,81,82,83,84,85], as well as their importance in catalysis [42,43,86].

In order to complement and confirm the conclusions obtained by the ETS-NOCV method we additionally plotted the contour of the reduced density gradient of the NCI (Non-Covalent Index) method [55] for KN(CH₃)₂BH₃, Figure 9. It was demonstrated that this method is well suited for visualization of non-covalent interactions [55]. It is clearly seen, in line with ETS-NOCV based study, that the crystal of KN(CH₃)₂BH₃ is stabilized by numerous inter and intramolecular K∙∙∙H–B as well as additionally by CH∙∙∙HC interactions. The same is valid for LiN(CH₃)₂BH₃. It is to be anticipated that the existence of strong M∙∙∙H–B interactions in both crystals might affect the mechanism of dehydrogenation similarly as has been already shown for the parent compound LiNH₂BH₃ for which dehydrogenation is initiated by Li∙∙∙H–B interactions and proceeds further through formation of LiH hydride and NH₂BH₂ as intermediates [87,88,89].

4. Conclusions

In the present study non-covalent interactions in the hydrogen storage materials LiN(CH₃)₂BH₃ and KN(CH₃)₂BH₃ are for the first time quantitatively (and qualitatively) described by means of the charge and energy decomposition method ETS-NOCV as well as the Interacting Quantum Atoms (IQA) approach.

It was found, in line with the pioneering work of McGrady et al. [29], that both crystals are stabilized by numerous intra- and intermolecular M∙∙∙H–B interactions (M = Li, K). The ETS-NOCV calculations indicated that these bonds are electrostatically dominated, followed by charge transfer and dispersion contributions. Interestingly, the intramolecular charge transfer contributing to Li∙∙∙H–B interaction appeared to be more pronounced than the corresponding intermolecular delocalization. We further noticed in LiN(CH₃)₂BH₃ that formation of intermolecular Li∙∙∙H–B interactions enforces charge delocalization within the homopolar BH∙∙∙HB contacts which explains the presence of the QTAIM bond critical points between clashing hydrogen atoms as was found by McGrady et al. [29]. In contrast, our investigations allowed us to conclude that monomers of LiN(CH₃)₂BH₃ are not likely to spontaneously form any stable aggregates via “pure” stabilizing BH∙∙∙HB interactions (BH₃ to BH₃ orientation) due to the presence of negatively charged borane units—from a fundamental point of view it demonstrates that BH∙∙∙HB interactions can be rather considered as destabilizing in this type of compounds. We further confirmed this conclusion by the Interacting Quantum Atoms (IQA) energy decomposition calculations—namely, the dihydrogen interaction energy in BH∙∙∙HB appared to be significantly positive (destabilizing), ${\mathrm{E}}_{\text{int}}^{\text{B}\left(\text{H}\right)\cdots \text{H}\left(\text{B}\right)}=+49.4\text{ kcal}/\text{mol}$ as opposed to Li∙∙∙H–B interactions, ${\mathrm{E}}_{\text{int}}^{\text{Li}\cdots \text{H}\left(\text{B}\right)}=-98.8\text{ kcal}/\text{mol}$.

Contrary to homopolar BH∙∙∙HB interactions, the ETS-NOCV and IQA methods allowed us to identify stabilizing homopolar dihydrogen interactions CH∙∙∙HC in both LiN(CH₃)₂BH₃ and KN(CH₃)₂BH₃—the presence of such stabilization has already been suggested by McGrady and coworkers based on a topological QTAIM study [29]. We found herein quantitatively that these interactions are dispersion dominated (64% for LiN(CH₃)₂BH₃ and 69% for KN(CH₃)₂BH₃), followed by charge transfer (13% for both LiN(CH₃)₂BH₃ and KN(CH₃)₂BH₃) and electrostatic (23% for LiN(CH₃)₂BH₃ and 17% for KN(CH₃)₂BH₃) terms. It was confirmed that these interactions are far stronger in the potassium crystal due to the larger number of CH∙∙∙HC contacts compared to the lithium analogue. Moreover, the NOCV-based deformation density contributions allowed to state that the “electronic” part of the CH∙∙∙HC interaction is based on both donation from the occupied σ(C–H) bonds into the empty σ*(C–H) of methyl groups as well as polarization of the C–H bonds (mixing of σ/σ*(C–H)).

Briefly summarizing, our in-depth theoretical investigations performed by means of the ETS-NOCV and IQA energy decomposition methods, electrostatic potentials and charges, allowed us to confirm most of the findings that have been already reported in the pioneering work of McGrady et al. [29]. Furthermore, we provided for the first time the energetic description of non-covalent interactions contributing to the stability of LiNMe₂BH₃ and KNMe₂BH₃ as well demonstrating, contrary to McGrady et al. [29], the repulsive nature of the homopolar interactions BH∙∙∙HB. The latter is in line with numerous experimental papers [68,69,70,71,72]. Due to the fact that our calculations are based on the cluster approach as well as the fact that all theoretical methods are not free from approximations and very often not from arbitrariness, we believe that further works are needed from both theoretical and experimental laboratories in order to fully uncover the nature of homopolar BH∙∙∙HB interactions.