s-Block metallabenzene: aromaticity and hydrogen adsorption

Second group metal dimers can replace the carbon atom in benzene to form metallabenzene (C5H6M2) compounds. These complexes possess some aromatic character and promising hydrogen adsorption properties. In this study, we investigated the aromatic character of these compounds using aromaticity indices and molecular orbital analysis. To determine the nature of interactions between hydrogen and the metallic center, variation-perturbational decomposition of interaction energy was applied together with ETS-NOCV analysis. The results obtained suggest that the aromatic character comes from three π orbitals located mainly on the C5H5 − fragment. The high hydrogen adsorption energy (up to 6.5 kcal mol−1) results from two types of interaction. In C5H6Be2, adsorption is controlled by interactions between the empty metal orbital and the σ orbital of the hydrogen molecule (Kubas interaction) together with corresponding back-donation interactions. Other C5H6M2 compounds adsorb H2 due to Kubas interactions enhanced by H2–π interactions. Graphical Abstract First π orbital in C5H6Be2 Electronic supplementary material The online version of this article (doi:10.1007/s00894-014-2552-6) contains supplementary material, which is available to authorized users.


Introduction
The concept of aromaticity was defined in the nineteenth century to explain the unusual properties of benzene.
Nowadays, this term is also used to describe other compounds including heterocyclic [1] and metallacyclic compounds [2][3][4] as well as metal clusters [5,6]. Since the term "aromaticity" is used for many different types of resonance stabilization phenomena, the classification of compounds as aromatic or nonaromatic is not always obvious. The quantitative determination of this property is possible via introduced aromaticity indices (AI). Several AI have been proposed; however, all of them are based on one of the following ring criteria: geometrical, energetic, magnetic or electronic properties. The geometrical aspect of aromaticity is measured predominantly by the HOMA (harmonic oscillator model of aromaticity) index proposed by Kruszewski and Krygowski [7]. HOMA is the most commonly used AI; however, this index depends on empirical parameters that are defined for organic compounds only. Therefore HOMA cannot be applied to organometallic compounds. NICS (nucleus-independent chemical shift) represents a series of indices [8] that allow evaluation of aromaticity on the basis of the magnetic properties of the ring. NICS uses a negative value of absolute shielding in the given position. NICS indices do not depend on empirical parameters and thus may be applied easily to any aromatic system. A negative value of an index indicates that the compound is aromatic whereas positive values states that the investigated compound is antiaromatic. There are three commonly used NICS indices: NICS(0), defined as negative value of NICS in the ring center; NICS(1), which represents negative value of NICS measured 1 Å above the center of a ring; and NICS(1)zz, which corresponds to anisotropic part of NICS (1) in the direction perpendicular to the ring. The ring center is usually defined as a geometric center; however, it can also be defined as a ring critical point (RCP) based on the atoms in molecules (AIM) theory introduced by Bader [9].
The key feature of aromatic compounds is a delocalization of π electrons, which leads to additional stabilization of a structure. The stabilization energy is called the resonance energy or aromatic stabilization energy (ASE). In general, ASE represents the energy difference between an aromatic compound and its non-aromatic unsaturated isomer. ASE can be calculated in the isodesmic reaction [10], where substrates and products have the same number of atoms, bonds, and bond types. The standard example constitutes the hypothetical reaction of three cyclohexenes that form benzene and two cyclohexanes. The other reaction suited to calculating ASE is a double-bond migration isomerization reaction [11]. In this reaction, the product constitutes an aromatic compound with a methyl substituent in the meta position whereas the substrate has a methylidene substituent and an additional hydrogen in the para position.
A number of indices are based on electronic properties. The aromaticity of a molecule can be characterized by properties of electron density calculated in RCP. Palusiak et al. [12] and Ebrahimi et al. [13] proved that the value of electron density, laplacian of electron density, and energies of electron density, including potential (V), kinetic (G) and total (H) electron energy, are strongly correlated with a HOMA value. Besides electronic properties in bond critical point (BCP), several other indices are defined within the AIM framework. However, most common indexes [14], e.g., PDI (para-delocalization index) [15], ΔDI [15], FLU [16] and FLU π [16], cannot be applied to C 5 M 2 H 6 compounds. The reasons are as follows [14]: by definition, the PDI index can be used only for 6membered rings; ΔDI requires a clearly defined Lewis structure; FLU, like HOMA, needs empirical parameters for each bond type whereas FLU π can be used only for planar compounds. One of few AIM-based indices that can be applied in the studied case is the Shannon aromaticity (SA) index [17] based on the Shannon entropy [18]. The formula describing SA is presented below: where with N representing the overall number of bonds in the ring and ρ BCP being electron density in the BCP. The value 0 corresponds to the homonuclear aromatic ring. SA increases when aromaticity decreases. Several metalloaromatic compounds have been reported in literature and the topic is covered by several reviews [2][3][4]. Most of these compounds are benzene analogues where CH groups are replaced by transition metals with ligands, so called metallabenzenes [19,20]. Compounds of this type containing transition metal atoms have been obtained experimentally and reveal aromatic properties (e.g., undergo aromatic substitution). Fernandez and Frenking [21] have proved that metallabenzenes containing Os, Ru, Ir, Rh, Pt, and Pd possess five π orbitals and can be considered as 10 electron Hückel aromatic systems.
Recently, we reported that the Be 2 fragment, which is isoelectronic with the carbon atom, can act as C(sp 2 ) [22]. The Be 2 center can replace the carbon atom in aromatic hydrocarbons and form compounds characterized by structure and molecular orbitals analogous to that in corresponding hydrocarbons. The aim of this study was to characterize the bonding scheme and aromaticity of compounds hosting Be 2 and other M 2 fragments (with M being second group metals; Table 1).
Compounds including a Be 2 center constitute promising materials for hydrogen storage, and can adsorb hydrogen with energy up to 6 kcal mol −1 . Hence, the second aim of this work was to characterize interactions of the M 2 moiety with hydrogen molecule(s). Here, we are concentrating on metallabenzene (C 5 M 2 H 6 ). Previous studies on Be 2 -containing compounds have shown that C 5 Be 2 H 6 preserves all features of higher aromatic hydrocarbons.

Computational details
If not mentioned otherwise, molecular structures and complexes with hydrogen molecule(s) were minimized at MP2 level of approximation [23] applying the aug-cc-pvdz atomic basis set [24]. Harmonic vibration frequencies revealed no imaginary frequencies indicating that structures correspond to true minima. NICS values were calculated using the gauge-independent atomic orbital (GIAO) method [25] applied at the same level of theory. Calculations mentioned above were performed in the Gaussian 09 software suite [26].
In order to investigate the nature of interactions between the M 2 center and molecular hydrogen, the total MP2 interaction energy was decomposed according to the hybrid variational-perturbational scheme [27]. In adopted scheme the total interaction energy is partitioned into Hartree-Fock (HF) (ΔE HF ) and second order perturbation (E MP2 ) terms: The HF interaction energy is further decomposed into first order electrostatic (E el ), Heitler-London exchange (E ex ) and higher order delocalization (E del ) terms. All energy contributions are calculated in the dimer-centered basis set and hence are free from basis set superposition error (BSSE). The interaction energy calculations were performed utilizing the GAMESS software [28] incorporating EDS modifications [29].
The interaction energy was also calculated applying the ETS-NOCV scheme [30] based on the ETS protocol [31]. In the ETS method, the total interaction energy (ΔE int ) consists of contributions from orbital interactions (E orb ), Pauli repulsion (E pauli ) and electrostatic interaction (E elst ): The NOCV approach [32] makes it possible to decompose deformation electron density into a set of NOCV densitydifferential pairs (charge flow channels) and to estimate the contribution of each pair to the orbital energy (E orb ). The ETS-NOCV decomposition was calculated based on B3LYP functional and ADZVP basis set [33] build from Slater functions using ADF software [34][35][36].
In both decomposition schemes, electrostatic and Pauli exchange terms describe the same interactions, namely, electrostatic attraction and Pauli repulsion between unperturbed electron density of monomers. E orb , E ex and E MP2 terms describe the energetic effect of electron density change between separate monomers and the complex. The ETS-NOCV scheme is implemented for DFT methods only and obtained energies are not corrected for basis-set superposition error. Therefore we use the hybrid variational-perturbational scheme for the decomposition and ETS-NOCV for further characterization of orbital interactions only. Energy values in both decomposition schemes are not equal due to the different level of theory and basis sets type (Gaussian vs Slater function) but are in close agreement (typically below 1 kcal mol −1 ), thus allowing the qualitative discussion of orbital interactions.
Contours of molecular orbitals (MO) were plotted using the Jmol program [37] whereas NOCV charge flow channels were visualized using ADF-GUI software [38]. Contour values were adjusted individually to provide optimal readability of the figures (values 0.05-0.07 a.u. where used for MOs, and 0.002-0.005 for charge flow channels). AIM analysis was performed using the MultiWfn program [39].

Results and discussion
"Be 2 -benzene" A beryllium dimer incorporated into carbon networks can be treated as a superatom. The Be 2 fragment is isoelectronic with the single carbon atom and its orbitals can mimic C(sp 2 ) hybridization [22]. This approach suggests that the Be 2 fragment should be involved in a π electron delocalization in C 5 Be 2 H 6 . On the other hand, "Be 2 -benzene" can be treated as a coordination compound with the C 5 H 5 − ligand. In such a situation, the π delocalization (and aromaticity) depends on the ligand and does not involve the Be 2 center. The second hypothesis is based on the work of Fernandez and Frenking [21], which states that [(C 5 H 5 )Rh(PH 3 ) 2 (Cl 2 )] can be treated as a coordination compound with C 5 H 5 − ligand. Molecular orbitals of "Be 2 -benzene" carring π electrons (Table 4) are similar to these in silabenzene and germabenzene (Table S3). The first MO possesses only a single nodal plane and spans over all the carbon atoms. The remaining two MOs possess two nodal planes; however, in contrast to C 5 H 6 Si and C 5 H 6 Ge, the last π MO in "Be 2 -benzene" shows significant density depletion on the heteroatomic center. The AIM analysis revealed (Table S5) that each beryllium atom is involved in bonding interactions with both neighboring carbons (denoted as C orto ) and nearest hydrogen atoms. The BCP was also found on the path connecting both C orto atoms. The laplacian in all BCPs is positive, indicating closed shell (ionic) interactions. The results of AIM analysis suggest that "Be 2 -benzene" has a partially ionic structure and can be considered as a coordinating compound with C 5 H 5 − and H − ligands around the Be 2 2+ center. Aromatic indices cannot give a clear answer about the aromaticity of "Be 2 -benzene". NICS(0) shows no aromaticity whereas NICS(1) (and especially NICS(1)zz) suggests significant aromatic character. Based on magnetic properties, the "Be 2 -benzene" molecule is less aromatic compared   to the C 5 H 5 − ligand alone (Table S1). Electron density properties in RCP and SA index state that "Be 2 -benzene" is the most aromatic compound of the C 5 H 6 M 2 series. The opposite conclusion can be drawn from ASE, which suggests that C 5 H 6 Be 2 is the least aromatic of the series. The Be 2 fragment can adsorb one hydrogen molecule with relatively high energy (5.1 kcal mol −1 ). The adsorption of first hydrogen is controlled by the Kubas-like interaction. The classical Kubas interaction [40] involves the σ-orbital of the hydrogen molecule and the empty d-orbital of the metal. Beryllium does not possess any d orbitals; instead, a porbital is involved in interactions. In the ETS-NOCV analysis, in addition to the Kubas-like interaction based on electron transfer from hydrogen to metal, a strong back-donation is observed. The adsorption of a second hydrogen on Be 2 due to lower energy of the Kubas-like interaction is relatively low (1.9 kcal mol −1 ).
The geometry (Table 2) and molecular orbitals (Tables 3, 4) of C 5 H 6 Mg 2 ("Mg 2 -benzene") are similar to those observed in "Be 2 -benzene". AIM analysis indicates the same bonding scheme. Hence, the compound can be also treated as a complex with an aromatic C 5 H 5 − ligand. The aromaticity of "Mg 2benzene" based on magnetic properties calculated 1 Å above the ring is comparable to the "Be 2 -benzene" case. In contrast to "Be 2 -benzene", NICS values for "Mg 2 -benzene" in RCP are negative, confirming the aromatic character of the compound. The electron density properties suggest lower aromaticity, whereas ASE suggests higher aromaticity in comparison to "Be 2 -benzene". The adsorption energy of first hydrogen molecule on "Mg 2 -benzene" is lower than that of "Be 2 -benzene" (3.8 kcal mol −1 ). The adsorption of this molecule is controlled by electrostatic and delocalization terms that constitute almost equal contribution and are two times higher than correlation effects. The geometry of the complex shows that the distances between magnesium and hydrogen atoms from adsorbed H 2 molecule are not equal. The first hydrogen atom (denoted as H1) is attached directly to Mg whereas the second (denoted as H2) is in contact with the Mg and C para atoms. This may suggest that the hydrogen molecule interacts via the empty orbital on magnesium and also with π-orbitals located on carbons. ETS-NOCV reveals that three deformation density channels are responsible for 94 % of the orbital energy ( Table 5). The first charge flow channel, which is responsible for 60 % of orbital interactions, indicates an increase in electron density on H1 and between magnesium and the H1 atom. The same charge flow channel shows the density depletion on H2 and on πorbitals of the carbon atom in the para and meta positions. The density change can be characterized as an electron transfer from the hydrogen molecule to an empty orbital on the metal and can be classified as a Kubas-like interaction. However, this interaction differs from a typical Kubas interaction and that observed in "Be 2 -benzene". The NOCV charge flow channels for "Mg 2 -benzene" shows that only one atom from H 2 is involved in the Kubas-type interaction. This effect is probably caused by hydrogen polarization. The electrons from the sigma-orbital are moved toward H1 and, as a consequence, the Kubas-like interaction involves only one hydrogen from the H 2 molecule. This interaction, which involves the p-orbital from the metal and one hydrogen atom will be referred to here as a "polarized Kubaslike" interaction. The second NOCV charge flow channel reveals density depletion on H1, and a density increase on H2 and π orbitals (mainly on carbon in the para position). The channel describes the interactions of hydrogen with π orbitals. The third channel shows an increase in density on the hydrogen molecule with very diffuse density depletion behind H 2 (not shown on the picture). This phenomenon can be characterized as a back-donation.
Adsorption of the second hydrogen on the Mg 2 center takes place on the second magnesium atom on the opposite site of the molecule but in a similar manner to the first. The adsorption energy is significantly lower compared to the adsorption of the first H 2 (2.3 kcal mol −1 ). The energy decomposition shows that absolute values of all energy components are lower and their proportions are changed slightly. The dispersion attractions are more important in binding of the second hydrogen, whereas the delocalization term is less significant. ETS-NOCV analysis showed the same charge flow channels but its contribution to the complex stabilization is lower. The adsorption energy of the second hydrogen on the same magnesium atom of "Mg 2benzene" is as low as 0.6 kcal mol −1 . The interaction energy decomposition reveals that the adsorption is controlled by electrostatic and delocalization terms. ETS-NOCV charge flow channels show that the delocalization term is dominated by the Kubas interaction and the corresponding back-donation.
In summary, the replacement of Be 2 by Mg 2 changes structure, orbitals, and aromaticity of compounds only slightly. However, the hydrogen adsorption mechanism on C 5 H 5 Mg 2 differs significantly from that on C 5 H 5 Be 2 .
"Ca 2 -benzene" The structure of "Ca 2 -benzene" (Table 2) differs significantly from that of "Be 2 -benzene" and "Mg 2 -benzene". The molecule no longer possesses a C 2 axis of symmetry and the calcium atoms are not symmetrically equivalent. The first calcium atom (denoted as Ca down ) is close to the plane and interacts with two neighboring carbon atoms, whereas the second calcium atom (Ca up ) lies above the C 5 H 5 plane. The Ca up atom also interacts with π orbitals of other carbon atoms. The AIM analysis (Table S5) shows close shell bonding interactions between Ca up and carbon in the para position (C para ). Due to this interaction, C para lies slightly above the C 5 H 5 plane. The orbital contours indicates that the bonding interactions between metal atoms and neighboring carbons and hydrogen are similar to these in "Mg 2 -benzene". However, the character of the third molecular orbital differs from that of "Mg 2 -benzene" and confirms that the π orbital(s) from other carbon atoms are involved in calcium binding. The calcium interactions with π orbitals destroy the symmetry and shape of the π orbitals. The first orbital spans all the carbon atoms; however, it has no planar symmetry. The second orbital covers C-C and C-H bonds and differs drastically from those in benzene. The third orbital possesses a localized character typical of isolated p orbitals. In summary, π orbitals possess partially delocalized character but the interaction with Ca up atom destroys their aromaticity.
Values of all NICS indices are close to zero indicating that the compound is not aromatic. In comparison to "Mg 2benzene", "Ca 2 -benzene" has lower electron density and lower laplacian of electron density in RCP, which confirms its lower aromatic character. However, the aromatic stabilization index leads to the opposite conclusion. The ASE value for "Ca 2 -benzene" is over two times higher than that for "Mg 2 -benzene". This leads to the interpretation that "Ca 2 -benzene" is more aromatic that "Mg 2 -benzene". The ASE index represents the difference between an aromatic compound and its non-aromatic isomer; normally, this energy corresponds to the aromatic stabilization. In the case of "Ca 2 -benzene", the calcium atom interacts with π orbitals, thus the ASE energy can be considered as the sum of aromatic stabilization and Ca-π interactions. Therefore, the picture drawn by NICS indices seems to be the most reliable. Hence, one can conclude that Ca 2 -benzene has negligible aromatic character.
Since calcium atoms in "Ca 2 -benzene" are not symmetrically equivalent, their interactions with hydrogen differ. The upper calcium adsorbs the first hydrogen with an energy of 3.8 kcal mol −1 , which is equal to adsorption of the first H 2 on "Mg 2 -benzene". The energy decomposition shows that the nature of adsorption is similar in both complexes; however, the dispersion component is less significant in the Ca 2 -benzene case. Two ETS-NOCV charge flow channels are responsible for 94 % of delocalization energy. The first channel shows a density increase on the hydrogen atom closer to Ca (H1) and between H1 and calcium. Electron depletion is observed on the hydrogen close to the aromatic ring (H2) and on the π orbitals of the carbon atom in the para position (C4) and in ortho positions. This channel also reveals the density increase between H2 and C4 atoms. The interactions constituting the first channel can be summarized as polarized Kubas-like, supported by H 2 -π interactions. The adsorption energy of the second H 2 on Ca up is 1.4 kcal mol −1 and electrostatic interaction constitutes the main attractive term. The charge flow channel reveals that the delocalization term (which is responsible for 36 % of attractive interactions) is controlled by the polarized Kubas interaction. The calcium, which is located close to the plane of the ring, can adsorb up to two hydrogen molecules but the energy is as low as 1.6 and 1.2 kcal mol −1 . The first H 2 is adsorbed above the aromatic ring and the interaction is controlled by electrostatic and dispersion contributions. The delocalization is responsible for only 22 % of attractive interactions. ETC-NOCV analysis indicated the lack of electronic interactions between hydrogen and the calcium atom. Such interactions were observed between H 2 and π orbitals. The adsorption of the second H 2 was controlled by electrostatic attractions, which were responsible for 49 % of attractive interactions, whereas delocalization and electron correlation contributions constituted 35 % and 16 %, respectively. The delocalization term is dominated by polarized Kubas, which is responsible for 89 % of the delocalization energy. "BeMg-benzene" The structure of "BeMg-benzene" is similar to that of "Be 2 -benzene" and "Mg 2 -benzene", i.e., all carbon and hydrogen atoms lie on a plane and the metal atoms are on the perpendicular axis. Distances of metal atoms from the plane are not equal: beryllium is located 0.5 Å below whereas magnesium lies 1.6 Å above it. In "BeMg-benzene", the Mg-C and Mg-H bonds are shorter compared to Mg 2 -benzene, which may indicate that the Mg 2 center is slightly too large to fit in the carbon atom position. However, replacement of Mg with the smaller Be atom does not affect the MOs. Orbitals involved in metal binding as well as π orbitals have the same character in "BeMg-benzene" as in pure compounds.
The NICS value in the RCP is close to zero, whereas 1 Å above RCP is negative, which is similar to the case for C 5 H 6 Be 2 and C 5 H 6 Mg 2 . Other AIs suggest that "BeMg-benzene" is more aromatic that "Mg 2 -benzene". Absolute values of electron density, density laplacian, and electron energy in RCP are approximately one-third higher that those of "Mg 2benzene". Similarly, the SA index of "BeMg-benzene" is onefourth lower compared to "Mg 2 -benzene" whereas ASE is over three times higher. Based on the results presented, one can conclude that the replacement of magnesium with beryllium in "Mg 2 -benzene" stabilized the structure and increased its aromatic character.
"BeMg-benzene" can adsorb the first hydrogen with relatively high energy of 6.5 kcal mol −1 . H 2 is located over the carbon ring close to magnesium. Adsorption is controlled by delocalization and electrostatic terms that are responsible for 44 % and 40 % of attractive interactions, respectively. The electron correlation constitutes only 16 %. The ETS-NOCV analysis revealed that three charge flow channels significantly stabilize the complex. The first, which can be characterized as polarized Kubas, is supported by interaction with π orbitals. This interaction is responsible for 60 % of delocalization energy. The second channel (22 % of delocalization energy) can be classified as a hydrogen interaction with π orbitals, whereas backdonation-the third channel-constitutes 13 % of overall delocalization energy. The adsorption energy of the second hydrogen on Mg amounts to 0.7 kcal mol −1 only. This interaction is also controlled by electrostatic and delocalization terms but details of the delocalization attraction are different. The main delocalization component (62 % of overall energy) is represented by Kubas interactions involving both hydrogen atoms (although not in equal proportions). The back-donation, which causes an increase of density on the adsorbed hydrogen, is responsible for 24 % of delocalization energy.
The hydrogen adsorption on beryllium in "BeMg-benzene" is low; the first H 2 can be bound with an energy of 1.2 kcal mol − 1 , whereas the second is bound with 0.2 kcal mol −1 only. Adsorption is controlled by electron correlation (46 % of attractive contribution) and electrostatics (37 %). Analysis of the delocalization term shows no interaction with the metal center, which rather comes from interactions with π orbitals. Similarly, the second hydrogen is bound by delocalization and electrostatic terms. It is worth noting that the hydrogen adsorption properties of beryllium in "BeMg-benzene" decreases dramatically compared to those of "Be 2 -benzene". A possible explanation for this fact is the high electron density (which can be seen on the bonding orbital) and low atomic charge on the Be atom. The NBO charge on Be in "BeMg-benzene" is only 0.540, whereas in "Be 2 -benzene" it is 0.822 electron.
"BeCa-benzene" and "MgCa-benzene" In "BeCa-benzene" and "MgCa-benzene" compounds, the smaller metal atom is close to the C 5 H 5 plane. Being above the ring, the calcium atom interacts not only with neighboring carbons but also with π electrons of other carbon atoms. The carbon-calcium bond length is 2.527 Å in "BeCa-benzene" and 2.507 Å in "MgCa-benzene". The interatomic distances between calcium and carbon in the meta position for "BeCabenzene" and "MgCa-benzene" are 2.705 Å and 2.725 Å, whereas the distances to carbon in the para position are 2.682 Å and 2.710 Å, respectively. AIM analysis in both cases indicates Ca-C para bonding interactions. Bonding orbitals in "BeCa-benzene" have no π character whereas the third orbital involved in bonding metals has a partial π character. The contour of π orbitals indicates that, in both compounds, the first orbital is delocalized over all carbons. The second π orbital in both compounds has one nodal plane, whereas the symmetry and shape of both compounds is different. "BeCabenzene" possesses a typical π orbital that spans symmetricaly over carbon atoms in ortho and meta positions, whereas in "MgCa-benzene" this orbital also covers the C-H bond. The third orbital in both compounds reveals interactions with metal. NICS indices show a lack of aromaticity in both compounds.
Like "Ca 2 -benzene" and "BeMg-benzene", the adsorption of hydrogen on "BeCa-benzene" and "MgCa-benzene" proceeds rather on the heavier metal located above the carbon ring. The calcium atom in "BeCa-benzene" can bind two hydrogen atoms with energies of 4.5 and 2.5 kcal mol −1 , respectively. Both interactions are controlled by electrostatic and delocalization terms. Orbital interactions reveal that 64 % of the delocalization contribution comes from polarized Kubas interactions supported by interactions with the π orbital of carbon in the para position, and 31 % from the "classical" Kubas orbital overlapping. The second hydrogen is located above the Ca-H bond. ETS-NOCV charge flow channels reveal that 93 % of the delocalization interaction came from polarized Kubas interactions. The mechanism of hydrogen adsorption on the Ca atom in "MgCa-benzene" is identical but the overall adsorption energy is lower by 3.9 and 1.9 kcal mol −1 , respectively.
The lighter element (Be and Mg) in the compounds presented can adsorb up to two H 2 molecules but with interaction energies as low as 1.1-1.4 kcal mol −1 . In "BeCa-benzene", the hydrogen adsorption energy is virtually the same for the first and second H 2 molecule; however, the energy components reveals a dramatically different mechanism. The energy decomposition (Table 6) indicates significant delocalization and electrostatic interactions (both above 10 kcal mol −1 ) between the first H 2 and "BeCa-benzene". ETS-NOCV charge flow channels  show that the delocalization interaction comes from Kubas interactions involving both hydrogen atoms (58 % of delocalization energy) and two polarized Kubas interaction between one hydrogen atom from the H 2 molecule and beryllium (22 % and 16 % of the delocalization term). Adsorption of the second hydrogen has a similar overall energetic effect (1.3 kcal mol −1 ) but the mechanism is different. The main attractive term comes from electron correlation, which is responsible for 48 % of the attractive energy, while the delocalization term constitutes only 17 % of overall attractions. ETS-NOCV decomposition shows no significant interactions between hydrogen and beryllium, and indicates that the delocalization term comes from interactions with π orbitals. A similar situation is observed in "MgCa-benzene", where adsorption of the first hydrogen is controlled by electrostatic and delocalization terms. The delocalization component comes from Kubas interactions supported by interactions with π orbitals (the first charge flow channel) and polarized Kubas interaction. Adsorption of the second Table 6 Interaction energy decomposition of hydrogen adsorption energy and the ETS-NOCV analysis.Definitions as in

Conclusions
Beryllium and magnesium can form metallabenzene-type compounds of general formula C 5 H 6 M 2 exhibiting some aromatic character. The aromatic character comes from the C 5 H 5 − fragment rather than from the whole molecule. The heavier group 2 metals (e.g., calcium) are too large to fit in the carbon position, hence C 5 H 6 Ca 2 has insignificant aromaticity. Aromatic properties of benzene analogues with two different metals depend on the atom size and relative difference in the size, hence the highest aromaticity is observed in C 5 H 6 BeMg, whereas compounds with calcium atom do not exhibit aromatic character.
Aromatic indices for such compounds should be interpreted with caution. The values obtained result not only from aromatic properties but also from the metal type. AIMbased AIs give preference to more electronegative atoms. The Table 7 The interaction of hydrogen(s) with "Be 2 -benzene". NOCV charge flow channels (blue/red contours represent accumulation/depletion electron density), corresponding energy (in kcal mol −1 ), and flow type. The flow channels are ordered by decreasing energy aasymmetric interaction, ssymmetric interaction, M←H 2electron transfer from hydrogen molecule to metal atom, MH←H 2electron transfer to MH bond, M←H 2 →Celectron transfer to metal and carbon atoms, C←H 2electron transfer to carbon atom(s), Back-donationelectron transfer to H 2 molecule ASE index for the compounds presented measure the sum of aromatic stabilization and metal-π interaction rather than the aromaticity alone. C 5 H 6 M 2 -type compounds can adsorb a single hydrogen molecule with energy up to 6.5 kcal mol −1 , depending on the metal type. Such a value is sufficient for commercial applications. Energy decomposition shows that adsorption is controlled by delocalization and electrostatic contributions. The adsorption mechanism in Be 2 -benzene involves only metal atoms; however, in other compounds the carbon atom also takes part in the hydrogen binding. NOCV charge flow channels reveal that one hydrogen atom from the H 2 molecule interacts with π orbitals of the C 5 H 5 − ligand whereas the second hydrogen atom interacts with an empty orbital of the metal. The aromatic character of compounds with two different metals (C5H6MM') lies between properties of systems with the same metal type (C 5 H 5 M 2 and C 5 H 5 M' 2 ). However, the interaction with hydrogen for such compounds is stronger than that for C 5 H 5 M 2 or C 5 H 5 M' 2 and adsorption occurs on heavier metals. This tendency is visible, especially in Table 8 The interaction of hydrogen(s) with "Mg 2 -benzene". NOCV charge flow channels (blue/red contours represent accumulation/depletion electron density), corresponding energy (in kcal mol −1 ), and flow type. The flow channels are ordered by decreasing energy aasymmetric interaction, ssymmetric interaction, M←H 2electron transfer from hydrogen molecule to metal atom, MH←H 2electron transfer to MH bond, M←H 2 →Celectron transfer to metal and carbon atoms, C←H 2electron transfer to carbon atom(s), Back-donationelectron transfer to H 2 molecule Table 9 The interaction of hydrogen(s) with "Ca 2 -benzene". NOCV charge flow channels (blue/red contours represent accumulation/depletion electron density), corresponding energy (in kcal mol −1 ), and flow type. The flow channels are ordered by decreasing energy aasymmetric interaction, ssymmetric interaction, M←H 2electron transfer from hydrogen molecule to metal atom, MH←H 2electron transfer to MH bond, M←H 2 →Celectron transfer to metal and carbon atoms, C←H 2electron transfer to carbon atom(s), Back-donationelectron transfer to H 2 molecule Table 10 The interaction of hydrogen(s) with "BeMg-benzene". NOCV charge flow channel (blue/red contours represent accumulation/depletion electron density), corresponding energy (in kcal mol −1 ), and flow type. The flow channels are ordered by decreasing energy aasymmetric interaction, ssymmetric interaction, M←H 2electron transfer from hydrogen molecule to metal atom, MH←H 2electron transfer to MH bond, M←H 2 →Celectron transfer to metal and carbon atoms, C←H 2electron transfer to carbon atom(s), Back-donationelectron transfer to H 2 molecule Table 11 The interaction of hydrogen(s) with "BeCa-benzene". NOCV charge flow channel (blue/red contours represent accumulation/depletion electron density), corresponding energy (in kcal mol −1 ) and flow type. The flow channels are ordered by decreasing energy aasymmetric interaction, ssymmetric interaction, M←H 2electron transfer from hydrogen molecule to metal atom, MH←H 2electron transfer to MH bond, M←H 2 →Celectron transfer to metal and carbon atoms, C←H 2electron transfer to carbon atom(s), Back-donationelectron transfer to H 2 molecule Table 12 The interaction of hydrogen(s) with "MgCa-benzene". NOCV charge flow channel (blue/red contours represent accumulation/depletion electron density), corresponding energy (in kcal mol −1 ) and flow type. The flow channels are ordered by decreasing energy aasymmetric interaction, ssymmetric interaction, M←H 2electron transfer from hydrogen molecule to metal atom, MH←H 2electron transfer to MH bond, M←H 2 →Celectron transfer to metal and carbon atoms, C←H 2electron transfer to carbon atom(s), Back-donationelectron transfer to H 2 molecule C 5 H 6 BeMg. In "BeCa-benzene" and "MgCa-benzene", this effect is rather low and the adsorption energy is comparable to that of "Ca 2 -benzene".