Study of Beryllium, Magnesium, and Spodium Bonds to Carbenes and Carbodiphosphoranes

The aim of this article is to present results of theoretical study on the properties of C⋯M bonds, where C is either a carbene or carbodiphosphorane carbon atom and M is an acidic center of MX2 (M = Be, Mg, Zn). Due to the rarity of theoretical data regarding the C⋯Zn bond (i.e., the zinc bond), the main focus is placed on comparing the characteristics of this interaction with C⋯Be (beryllium bond) and C⋯Mg (magnesium bond). For this purpose, theoretical studies (ωB97X-D/6-311++G(2df,2p)) have been performed for a large group of dimers formed by MX2 (X = H, F, Cl, Br, Me) and either a carbene ((NH2)2C, imidazol-2-ylidene, imidazolidin-2-ylidene, tetrahydropyrymid-2-ylidene, cyclopropenylidene) or carbodiphosphorane ((PH3)2C, (NH3)2C) molecule. The investigated dimers are characterized by a very strong charge transfer effect from either the carbene or carbodiphosphorane molecule to the MX2 one. This may even be over six times as strong as in the water dimer. According to the QTAIM and NCI method, the zinc bond is not very different than the beryllium bond, with both featuring a significant covalent contribution. However, the zinc bond should be definitely stronger if delocalization index is considered.

Although Hund's rule favors the spin triplet state over the singlet one [146], the requirements that invert this relationship, i.e., make the singlet state an electronic ground state, are well known. This may happen if either some appropriate geometric requirements are met [147][148][149] or one or both of the substituents R 1 and R 2 are σ-electron-withdrawing or π-electron-donating [149][150][151][152][153]. The latter requirement is met especially in the presence of strongly electronegative atoms with lone electron pairs, such as P, N, O, F, Cl, etc. In this case, the preference for the singlet state results from partial delocalization of the electron charge from electron lone pairs of these atoms to the unfilled p orbital on the carbene carbon atom (Figure 1).
Apart from carbenes, an equally important and interesting group of organic compounds is the so-called carbodiphosphoranes (CDPs) and their amine analogues [154][155][156][157][158][159][160][161][162][163][164][165][166][167][168][169][170][171][172]. Their uniqueness in the electronic structure (see Figure 1) is that, unlike the previously described carbenes, in CDPs, none of the four valence electrons of the carbon atom participate in ligand binding, and therefore these electrons remain unbound. Instead of covalent bonds as in carbenes, the carbon atom in CDPs is bound to ligands via donor-acceptor R→C bonds [167]. These non-binding valence carbon electrons form two lone pairs, and not just one as in singlet carbenes. It should therefore be expected that CDPs exhibit greater nucleophilic abilities than singlet carbenes, and moreover, they should be felt not only in the plane of the molecule but also in the direction perpendicular to it.
It is understandable that so far, the vast majority of theoretical studies on beryllium and magnesium bonds have used as Lewis bases small molecules containing either some atoms with good electron-donating properties or π bonds [23][24][25][26][27][29][30][31][32]. On the other hand, reports of systems containing spodium bonds [91][92][93], especially with zinc (they could be called zinc bonds) are very rare [91,93]. It is also quite understandable that the research on carbenes and CDPs is mostly experimental. This is of course due to their huge role in organic and organometallic synthesis, as mentioned earlier. For this reason, beryllium bonds, magnesium bonds, or spodium bonds with the zinc atom as the Lewis acid center (i.e., the zinc bonds) with the participation of either carbenes or CDPs are most often found by crystallographic methods in the solid state. In this case, both the carbene (or the CDP) and the Lewis acid interacting with it are molecules containing many different substituents and functional groups, often of considerable size, which makes the systems themselves also generally bulky.
In order to unite these two thematic areas, this article describes the result of theoretical research on a large group of dimers with a beryllium bond, magnesium bond, or zinc bond between various Lewis acids of the MX 2 (where M = Be, Mg, Zn and X = H, F, Cl, Br, Me) type and some fundamental carbenes ((NH 2 ) 2 C, imidazol-2-ylidene, imidazolidin-2ylidene, tetrahydropyrymid-2-ylidene, and cyclopropenylidene) and CDPs ((PH 3 ) 2 C and (NH 3 ) 2 C) acting as a Lewis base. Therefore, the aim of this article is to present the results of theoretical research on the properties of C· · · M bonds, where C is either a carbene or CDP carbon atom. It should be noted that due to the aforementioned scarcity of reports, in particular theoretical ones, on systems featuring a zinc bond (i.e., the spodium bond [1] with the participation of a zinc atom acting as a Lewis acid center), the reported studies on the C· · · Zn bond-containing systems investigated here represent an especially considerable novelty. At the same time, the presented results on the properties of this bond and slightly similar C· · · Be and C· · · Mg bonds will contribute to increasing the knowledge of both the carbenes chemistry and the chemistry of CDPs. It is worth mentioning at this point that the presence of X halogen atoms leads in some of the dimers considered here to certain symptoms that indicate interactions accompanying the leading C· · · M bond. Therefore, one of the sub-goals of this article is to investigate the conditions that favor these additional weak interactions.

Results and Discussion
As mentioned in the Introduction, this article describes research on systems containing a beryllium bond, a magnesium bond or a zinc bond, where the role of Lewis acids is played by the MX 2 molecules (where M = Be, Mg, Zn and X = H, F, Cl, Be, Me), while the role of Lewis base is played by either carbene ((NH 2 ) 2 C, imidazol-2-ylidene, imidazolidin-2-ylidene, tetrahydropyrymid-2-ylidene, or cyclopropenylidene) or carbodiphosphorane (being either (PH 3 ) 2 C or (NH 3 ) 2 C). The monomers themselves and their dimers, in which the described bonds occur, are presented in separate subsections.

Investigated Systems
The considered MX 2 molecules are characterized by a linear structure in which M-X bonds are formed by overlapping of the hybridized sp orbital of the metal atom with one of the orbitals of X. Due to lower electronegativity of the metal atom, this atom is endowed with a partial positive charge (Table 1), becoming electron-depleted and therefore a Lewis acid center. The atomic charge values shown in Table 1 confirm the known fact that they can be very significantly dependent on the method of obtaining them in the calculations [173][174][175][176]. The atomic charges obtained by the NBO and QTAIM methods seem to be greatly exaggerated. In the context of the presented results, however, it is more important that all the methods of obtaining atomic charges used here (i.e., Hirshfeld [177][178][179], NBO [180,181], and QTAIM [182][183][184]) show that in the set of MX 2 molecules for a given metal M the most positive charge on the M atom occurs when X = F. This is fully understandable due to the very high electronegativity of the fluorine atom. Conversely, the smallest positive charge on the M atom occurs when X = H. This result is not as expected, because, due to the positive inductive effect (+I) of the methyl group, one would expect the smallest positive atomic charge of M in MMe 2 . It is also seen that the Cl and Br atoms lead to similar atomic charges on M. Importantly, all the methods used show that the highest positive charge occurs in MgF 2 , and the lowest in ZnH 2 . If we refer to the most reliable [175,176] Hirshfeld atomic charges, then these values are 0.924 and 0.354 au, respectively. The former value suggests an extremely high polarization of the Mg-F bond, which practically becomes the Mg + F − ionic one. A practical consequence of this finding is that, assuming electrostatic reasoning, the MgF 2 molecule should be the best Lewis acid, and therefore it should theoretically form the strongest adducts with carbenes and CDPs.
With the values of the atomic charges obtained by various theoretical methods, it is interesting to see if there are clear relationships between them. Figure 2 shows the relationships between the Hirshfeld charges and their equivalents obtained by the NBO or QTAIM method.  Table 1) obtained by the Hirshfeld method and either NBO or QTAIM.
As can be clearly seen, the linear relationships between the Hirshfeld atomic charges and those obtained by the NBO or QTAIM method are very weak. Particularly in the case of the latter method, the obtained coefficient of determination is unacceptably low. This result shows that especially the atomic charges obtained by QTAIM should not be treated as reliable. This flaw of QTAIM-based atomic charges was attributed to irregular shapes of atomic basins, which give them multipolar moments.
The electrophilic properties of a particular metal atom, which is an acidic center in the MX 2 molecule, can be nicely illustrated by means of maps of the distribution of the molecular electrostatic potential (MESP) projected onto the electron density isosurface, as shown in Figure 3. The use of same scale of the electrostatic potential values (from 0.0 au (red) to 0.2 (blue)) for all MX 2 molecules allows one to easily capture the existing relationships. It is clearly seen that, upon going in the series Me→H→Br→Cl→F, i.e., from left to right in Figure 3, a belt of even more positive electrostatic potential develops around the central metal atom. This is of course confirmed by the corresponding values of the maximum electrostatic potential on M (V max (M)), which are provided in the last column of Table 1. For zinc compounds, these values (in au) increase in this series as follows: 0.037 < 0.050 < 0.065 < 0.071 < 0.093. Although the MESP maps for zinc molecules are very similar to those for beryllium, it is worth noting that in the former case, the corresponding MESP belt is wider and larger in diameter due to the larger atomic radius of Zn 2+ (88 pm) compared to Be 2+ (59 pm) [185]. Therefore, compared to beryllium, the zinc atom should be more accessible. The more important result, however, is that, for a given X, the belts of positive MESP are most visible when the central metal atom is magnesium. The values of V max (M) increase monotonically quite quickly in the order given earlier, reaching a maximum value of 0.243 au in MgF 2 ( Table 1). The fact that the V max (M) values increase in this order, while q(M) does not, suggests that V max (M) is perhaps the more reliable parameter describing the acidic nature of the central metal atom in MX 2 molecules than q(M). Although the linear relationship between the value of V max (M) is not very good (R 2 = 0.795) either when q(M) is computed utilizing the Hirshfeld method, it is much better than in the case of NBO-and especially QTAIM-based charges ( Figure 4).

• Carbenes and CDPs
Imidazol-2-ylidene, imidazolidin-2-ylidene, tetrahydropyrymid-2-ylidene, cyclopropenylidene, and (NH 2 ) 2 C have been used as model representatives of carbenes. In particular, the first two carbenes are often used in organic and organometallic chemistry and represent an important starting point in the syntheses of larger carbene compounds [98,102,103]. The CDPs group is represented by (PH 3 ) 2 C and its amino derivative (NH 3 ) 2 C. Both are the starting molecules for more complex CDPs obtained by substituting hydrogen atoms in one or both of the -PH 3 or -NH 3 groups. It is worth mentioning here that the fully saturated phenyl derivative, i.e., (PPh 3 ) 2 C was the first synthesized CDP [154]. Some fundamental parameters characterizing the considered carbenes and CDPs are presented in Table 2. Table 2. Some fundamental data characterizing carbenes and CDPs: the R-C-R angle (α RCR ), the atomic charge (Hirshfeld-, NBO-or QTAIM-based) of C (q(C) in au), the minimum value of the electrostatic potential on C (V min (C) in au), the energy of HOMO (E HOMO in eV)). When analyzing the obtained values of the atomic charge on the carbon atom, one can easily notice their great diversity, even in terms of sign. In the case of carbenes, positive QTAIM atomic charges have been obtained. Additionally, this method has given (too) large variation in the negative values on the C atom in (PH 3 ) 2 C and (NH 3 ) 2 C (−2.261 and −0.179 au, respectively). Both of these findings strongly suggest that atomic charges of QTAIM are highly unreliable. A similar conclusion applies to the atomic charges of NBO, although the values themselves are not that large. It is worth mentioning that the value of the atomic charge on a carbon atom of −1.43 au in (PPh 3 ) 2 C was used by Tonner et al. [167] as an argument supporting the bonding scheme of CDPs presented in Figure 1. However, taking into account large dependence of the atomic charge on the method used in calculations, it seems that this argument was perhaps not entirely correct. The more reliable [175,176] Hirshfeld atomic charges are negative in both carbenes and CDPs. Understandably, in the latter case they are much larger, which results from the role of the carbon atom as an acceptor in the R→C bonds ( Figure 1).

Molecule
Further valuable information on the nucleophilic abilities of singlet carbenes and CDPs can be obtained from the values of the minimum electrostatic potential on C (the penultimate column in Table 2) and the distribution of this potential around this atom (see Figure 5). The electrostatic potential maps clearly show the negative potential area around the C(2) atom in the carbenes or the C(0) atom in the CDPs. On the other hand, strong positive potential concerns mainly hydrogen atoms in strongly polar N-H bonds. While the characteristics of the negative potential distribution around the carbon atom are similar in carbenes (which is in line with the rather similar values of V min (C); Table 2), there is a clear difference between (PH 3 ) 2 C and (NH 3 ) 2 C. Specifically, in the latter case, this area in much clearer and much more spread around the carbon atom, which better emphasizes the great nucleophilic properties of this molecule. Both of these molecules also differ considerably in the value of V min (C) (−0.067 and −0.109 au, respectively).
Further information on the reactivity of molecules can be obtained from the Frontier Molecular Orbital theory [186,187], which has found its mathematical support in the Klopman-Salem Equation [188,189]. According to it, the electron-donating properties of the molecule can be characterized by the energy of HOMO. These energies for carbenes and CDPs are shown in the last column of Table 2. By far the least negative value of the HOMO energy obtained for (NH 3 ) 2 C (−4.35 eV) confirms that this molecule should undoubtedly be the most reactive, willingly acting as a Lewis base. It should be noted, however, that the HOMO energy, like the LUMO energy, which is also often used in the Frontier Molecular Orbital theory, is a global quantity, i.e., resulting from the electronic structure of the entire molecule, and therefore it does not necessarily correctly assess the nucleophilic and electrophilic properties of a molecule, which are most often strongly local. Moreover, these energies do not necessarily correlate well with the parameters characterizing the dimer strength. For example, as shown by Martín-Sómer et al. [24], LUMO energies correlate well with interaction energies (of some beryllium bonds) only when they are computed for acceptor molecules in their dimer geometries. For this reason, LUMO energy values for the fully optimized MX 2 molecules were not exposed in Table 2. Moreover, in the case of MX 2 molecules, the LUMO energy strongly depends on the X-M-X angle (α XMX ), decreasing considerably with increasing deviation from the linearity of the molecule. In this way, Martín-Sómer et al. [24] explained the large non-linearity of the BeH 2−n X n (X = F, Cl, Br; n ≤ 2) molecules in their dimers with ammonia. Therefore, it seems that there is nothing to prevent the same cause of MX 2 bending also working for other Lewis bases, such as the carbenes and CDPs considered here. It is also worth mentioning that the electron lone pair, which in carbenes is HOMO (quite strongly delocalized), in the case of CDPs, i.e., (PH 3 ) 2 C and (NH 3 ) 2 C, becomes HOMO-1, while HOMO is the electron lone pair perpendicular to the plane of the molecule ( Figure 6).

Dimers
The previous subsection has shown that in MX 2 (M = Be, Mg, Zn; X = H, F, Cl, Br, Me) molecules, the metal atom is a relatively strong acid center, while the C(2) atoms in the carbenes and C(0) in the CDPs are strong basic regions. Moreover, these atoms are the only such regions in these molecules (see Figures 3 and 5). Due to this alignment in electronic properties, it should be expected that the MX 2 molecules quite easily form a M· · · C bond to the C(2) carbon in carbenes or C(0) in CDPs. If so, it should lead to a particular type of beryllium, magnesium, or zinc (spodium) bond. As mentioned in the Introduction, the main purpose of this article is to describe these interactions. Nevertheless, the electrostatic potential distributions for MX 2 ( Figure 3) and carbenes and CDPs ( Figure 5) suggest that other interactions accompanying the leading M· · · C interaction may also be possible. In particular, some symptoms of the presence of a hydrogen bond of the N-H· · · X type (where X is a halogen atom, especially F) are to be expected. The geometries of the fully optimized dimers are shown in Figure 7. It is convenient to describe the characteristics of the systems containing carbenes and CDPs separately.

• Carbene dimers
The basic parameters characterizing the investigated carbene-containing dimers are shown in Table 3.  Table 3. Some fundamental data characterizing carbene· · · MX 2 dimers: C· · · M distance (in Å), changes of MX1 and MX2 bond lengths (in Å), XMX, CMX1, CMX2, LCL angles (in degrees), dissociation energy (in kcal/mol), charge transfer (in au).  Due to the simple structure of the cyclopropenylidene molecule, dimers containing this carbene will be discussed first. It should be noted that the plane of the slightly bent MX 2 molecule is oriented perpendicular to the plane of the cyclopropenylidene ring ( Figure 7). For this reason, the interaction between MX 2 and cyclopropenylidene is free from any other significant interactions than C· · · M. Although the earlier analysis of the values of atomic charges and electrostatic potentials on M and C suggests that the strongest C· · · M interaction should be present in the case of MgF 2 and the weakest in the case of ZnH 2 , this is not in line with the values of the distance C· · · M (d C···M ). Rather, these distances result from the radius of the metal atom, so in the case of beryllium, d C···M is less than ca. 1.83 Å, while in the case of Mg and Zn, this distance is over 2 Å. The penultimate column in Table 3 shows that cyclopropenylidene· · · MX 2 dimers are strongest (32-35 kcal/mol) when M is either Be or Mg and X is a halogen, especially Cl or Br. The lowest dissociation energy (10.3 kcal/mol) has been obtained in the case of ZnMe 2 . The weakest C· · · M in the presence of methyl groups has also been obtained in the case of M = Be or Mg and is in line with the weak +I character of the methyl group. Due to the C 2V symmetry, the following relations hold: ∆d MX1 = ∆d MX2 = ∆d av MX and α CMX1 = α CMX2 . The greatest elongation of the MX bond (0.093 Å) occurs in BeBr 2 . Along with a similar BeCl 2 , in this molecule, there is also the greatest deviation from linearity (α XMX amounts to ca. 134 • ). Thus, the geometric characteristics of the MX 2 molecule itself and the obtained D 0 values suggest that in the cyclopropenylidene· · · MX 2 dimers, the interaction should be strongest for BeBr 2 and BeCl 2 . It is interesting to see if similar finding also apply to dimers involving the other carbenes.
As Figure 7 shows, the MX 2 molecule lies in the same plane as the backbone atoms of the carbene molecule. This arrangement is also characteristic for dimers involving CDPs, (PH 3 ) 2 C and (NH 3 ) 2 C. In at least some cases, the planar geometry of the dimer can be explained by additional beneficial interactions (as will be discussed). As was the case with cyclopropenylidene, the intermolecular distance C· · · M is much shorter for beryllium (ca. 1.76-1.85 Å) than for either magnesium (ca. 2.17-2.28 Å) or zinc (2.03-2.20 Å). However, this does not mean stronger C· · · M interactions. The calculated dissociation energy values clearly show ( Table 3) that, as was the case for cyclopropenylidene, the strongest intermolecular C· · · M interaction occurs for BeBr 2 and BeCl 2 . Although the bond strength of the former is ca. 47-48 kcal/mol, it reaches up to 53 kcal/mol when BeBr 2 interacts with tetrahydropyrymid-2-ylidene. On the other hand, similar to cyclopropenylidene, the C· · · M interaction is the weakest (but clearly stronger than that of cyclopropenylidene) when the MX 2 molecule is ZnMe 2 . Consequently, in the dimers considered here, the dissociation energies of C(2)· · · M have a wide range from 10 to 53 kcal/mol. This result is in full accord with the recent generalization given by Alkorta and Legon that beryllium and magnesium bonds (the current results show also include the zinc bonds) are generally much stronger than hydrogen bonds, halogen bonds, etc. [29].
The LCL angle change (∆α LCL ) values show that the interaction between MX 2 and the carbene molecule leads to the opening of the latter molecule, with the effect being the greatest for (NH 2 ) 2 C (e.g., ∆α LCL = 6.8 • for (NH 2 ) 2 C· · · BeCl 2 ). This shows that the α NCN angle in (NH 2 ) 2 C is more flexible than in cyclic and therefore more rigid imidazol-2-ylidene, imidazolidin-2-ylidene and tetrahydropyrymid-2-ylidene (∆α CCC in cyclopropenylidene is negligible). Although in general the α NCN angle-opening effect in the carbene molecule does not seem to be dependent on the strength of the interaction with MX 2 , such a relationship can be found when comparing systems with similar skeleton stiffness. Therefore, in the group of the aforementioned cyclic carbenes, the strongest effect occurs in tetrahydropyrymid-2-ylidene (4.8 • ). Excellent linear relationships have been found (see Figure 8) between the change of the opening angle α LCL and the dissociation energy of the carbene· · · MX 2 dimer as long as the carbenes and the MX 2 molecules are treated separately. Note that the greater sensitivity of the opening angle in the case of the (NH 2 ) 2 C carbene is evident here by slightly larger slopes of the corresponding (red) lines. Moreover, the slopes of the linear fits for cyclic carbenes are similar to each other.
A characteristic effect that occurs during an interaction of the initially linear MX 2 molecule with a strong Lewis base is its significant bend [23]. For example, Martín-Sómer et al. have reported XCX angles (α XMX ) of 134 • -139 • , (B3LYP/6-311+G(3df,2p)) for dimers of X-substituted (X = F, Cl, Br) BeX 2 derivatives with ammonia [24]. This bending effect is much less (138 • -149 • ,) in BeX 2 (X = H, F, Cl) dimers with ethylene or acetylene, being much weaker Lewis bases interacting via π bonds [26]. The high sensitivity of the α XMX angle makes it particularly interesting to trace its values in the considered dimers. Due to the large number of the studied set of systems and their diversity (different acid centers M, different X substituents, different carbenes), a fairly wide range of α XMX variability has been obtained, from 131 to 151 • , i.e., as much as 20 • . The bending effect is greatest for BeCl 2 and BeBr 2 and the smallest for ZnMe 2 . The linear correlation between the XMX angle and the dissociation energy is acceptable for ZnX 2 (Figure 9, left) and the dimers of either imidazol-2-ylidene (R 2 = 0.942) or imidazolidin-2-ylidene (R 2 = 0.922) with BeX 2 (not shown). The fitting line for cyclopropenylidene has slightly different slope than the other four cases, which may result from different (perpendicular) orientation of the interacting molecules (Figure 7). The weak linear correlation for the remaining cases of carbene· · · MX 2 (M = Be, Mg) dimers may, at least partly, result from the presence of additional interactions is some of the considered dimers, which should have some influence on the angle XMX. In the case of the dimers involving ZnX 2 , as a consequence of good linear relationships between ∆α LCL and D 0 ( Figure 8) and α XZnX and D 0 (Figure 9, left), one also observes good linear relationships between ∆α LCL and α XZnZ (Figure 9, right).    . Relationships between either the XZnX angle (α XZnX ) and the dissociation energy (D 0 ) (left) or the change of the LCL angle (∆α LCL ) and α XZnX (right) for the L 2 C· · · ZnX 2 dimers (where L 2 C = (NH 2 ) 2 C (red), imidazol-2-ylidene (green), imidazolidin-2-ylidene (blue), tetrahydropyrymid-2-ylidene (magenta), cyclopropenylidene (black)).
Another effect observed during the formation of the carbene· · · MX 2 dimers is a significant elongation of both MX bonds. It should be clearly underlined here that, in general, both MX bond elongations are not necessarily of equal magnitude, so it is not necessarily true that ∆d MX1 = ∆d MX2 = ∆d av MX (Table 3). These unsymmetrical elongations of MX result from the presence of certain accompanying interactions in some of the dimers studied here.
Such cases are also clearly visible from different values of CMX1 and CMX2 angles (α CMX1 and α CMX2 , respectively) in Table 3. In such cases, the smaller of these angles (α CMX1 ) takes a value roughly about the right angle.
As already mentioned, any significant additional interactions are impossible in cyclopropenylidene dimers. In this case, the effect of MX bond elongation is therefore symmetrical, which allows for straigtforward analysis of the obtained relationships. The greatest elongation of the MX bonds is for X = Cl or Br, but only when the M atom is either beryllium or zinc (up to 0.093 Å for BeBr 2 ). Hence, the elongation effect is not entirely consistent with the strength of C· · · M if measured by D 0 . On the contrary, the smallest elongation of the MX bond occurs in MgMe 2 (0.030 Å) and MgF 2 (0.031 Å). Although the relatively small magnitude of the effect in the former case can be explained by a relatively weak interaction (the largest distance C· · · M amounting to 2.288 Å and the smallest bending of 149.1 • ), BeH 2 is also characterized by a small MX bond elongation (0.039 Å), and this molecule forms the shortest contact with cyclopropenylidene, amounting to only 1.743 Å. For the latter molecule, i.e., MgF 2 , the small effect of the MgF bond elongation can most likely be explained by a high polarity of the bond and therefore its considerable resistance to changes. It seems that the magnitude of the MX bond elongation does not clearly depend on d C···M or the interaction strength measured by D 0 .
As mentioned earlier, the asymmetry of the MX elongation effects in case of many dimers involving the remaining carbenes makes the analysis much more difficult, but mean value (∆d av MX ) provides some information. Regardless of carbene, this value for BeBr 2 is always greater than 0.106 Å and reaches a maximum value of 0.111 Å for tetrahydropyrymid-2-ylidene, thus confirming that presumably the C· · · M interaction is the strongest in the tetrahydropyrymid-2-ylidene· · · BeBr 2 dimer. The occurrence of the minimum values of ∆d av MX appears to be more irregular. Although BeH 2 is generally characterized by low values (ca. 0.048 Å), the lowest values (ca. 0.040 Å) are nevertheless found for MgMe 2 interacting with either imidazol-2-ylidene or imidazolidin-2-ylidene.
The last column in Table 3 shows values of charge transfer calculated by means of the most reliable [175,176] Hirshfeld atomic charges (CT H ). First, it should be noted that the obtained values are negative, which means that the formation of the carbene· · · MX 2 dimer leads to an increase in the total charge on the MX 2 molecule. Secondly, the obtained values are very large. Suffice it to mention that the corresponding charge transfer values obtained (on the same level of theory) for dimers HOH· · · OH 2 and HOH· · · NH 3 are −0.098 and −0.122 au, respectively. Thus, even the weakest charge transfers obtained for the investigated dimers are over two times greater (e.g., −0.270 au for cyclopropenylidene· · · ZnMe 2 ) and even reach almost four times higher values in some dimers with BeBr 2 (e.g., CT H amounts to ca. −0.46 au for imidazol-2-ylidene, imidazolidin-2-ylidene, and tetrahydropyrymid-2-ylidene). Undoubtedly, therefore, the carbene· · · MX 2 (M = Be, Mg, Zn) dimers considered here are characterized by a significant charge transfer, which is particularly high in the presence of highly polarizable halogen atoms in MX 2 , especially Br. This finding is also manifested by very good (R 2 = 0.955) linear correlation between CT H and ∆d av.
MX when X = Br and only slightly worse for X = Cl (R 2 = 0.917), while this correlation is very weak (R 2 = 0.154) for much less polarizable fluorine ( Figure 10). For the cyclopropenylidene· · · MX 2 dimers, there are also very good linear correlations between CT H and either ∆α LCL or D 0 (in the latter relationship, except in the case of M = Be) if only systems with different M atoms are treated separately (Figure 11). Unfortunately, similar relationships are generally much worse for other carbenes, which can be explained by the presence of additional intermolecular interactions in some of them, which to some extent affects the obtained values of the analyzed parameters.

• CDPs dimers
The fundamental data characterizing CDP· · · MX 2 dimers are included in Table 4. Its penultimate column shows that the C(0)· · · M interactions in the dimers formed by (PH 3 ) 2 C are comparable in strength to the C(2)· · · M bonds formed by the investigated carbenes, whereas those formed by (NH 3 ) 2 C are much stronger. Again, the maximum value is found for BeBr 2 , reaching 84 kcal/mol, a value comparable to the energy of weaker covalent bonds [190]. The other dimers with high values of D 0 are (NH 3 ) 2 C· · · BeCl 2 and (NH 3 ) 2 C· · · ZnF 2 (ca. 79 kcal/mol). It is worth recalling here the theoretical research by Jabłoński and Palusiak [108] on the ability of carbenes and CDPs to form hydrogen bonds. The results of those studies have shown that for the same Lewis acid (e.g., HCCH), the hydrogen bond to (NH 3 ) 2 C is much stronger than to (PH 3 ) 2 C (the MP2/augcc-pVTZ-based BSSE-corrected interaction energies amount to −9.16 and −5.31 kcal/mol, respectively, [108]), which further confirms the greater basicity of the former molecule. Although the (NH 3 ) 2 C· · · BeBr 2 and (NH 3 ) 2 C· · · BeCl 2 dimers are characterized by short C· · · Be distances (1.643 and 1.661 Å, respectively), the short C· · · Be distance (1.655 Å) is also present in the (NH 3 ) 2 C· · · BeH 2 dimer with much weaker interaction (64 kcal/mol).
Very high bond strength of C(0)· · · M in the (NH 3 ) 2 C· · · MX 2 dimers is in line with high values of charge transfer, which can even reach -0.610 au in the (NH 3 ) 2 C· · · BeBr 2 dimer. This value is more than six times greater than that of the water dimer and exactly five times greater than that of the water-ammonia dimer. A curiosity is the relatively low CT H value (−0.287 au) obtained for the (NH 3 ) 2 C· · · ZnF 2 dimer with a simultaneous very high dissociation energy (78.8 kcal/mol). In the next subsection, however, it will be shown that this dimer is characterized by a highly advanced proton transfer from N to F, which results in the formation of the N· · · H-F hydrogen bond. The formation of the H-F bond requires some removal of the electron charge from the fluorine atom.
An interesting result is that, as in the case of carbenes (Table 3), the interaction between MX 2 and (NH 3 ) 2 C causes a significant opening of the α NCN angle, whereas in the case of (PH 3 ) 2 C, the change in the α PNP angle is much smaller and may have different sign, most often being negative. This finding clearly differentiates the nitrogen atom from the phosphorus atom.
The comparison of the values of ∆d MX1 and ∆d MX2 as well as α CMX1 and α CMX2 shows a clear difference between the dimers with (PH 3 ) 2 C and the dimers with (NH 3 ) 2 C. Specifically, the former of them are characterized by the equality of both quantities, which indicates symmetry of these dimers with respect to the axis passing through the C and M atoms (see also Figure 7). In the latter case, however, this symmetry is clearly broken in most of the dimers, which results from the presence of other interactions accompanying the leading contact C· · · M. As a result, the search for linear correlations between the parameters from Table 4 for systems with (NH 3 ) 2 C is pointless, while the search for such correlations for systems with (PH 3 ) 2 C seems to be justified. Indeed, some reasonable linear correlations have been found, such as, for example, between d C···M and D 0 (see Figure 12) when the metal atom is either Mg (R 2 = 0.956) or especially Zn (R 2 = 0.989). On the other hand, when the acidic metal center is beryllium, the linear correlation clearly deteriorates (R 2 = 0.795). At least in part, this may be due to the much shorter Be-X bond compared to the Mg-X or Zn-X bonds, and thus stronger, although still rather weak, intermolecular interactions of the type -(PH 3 )· · · X. Moreover, quite good linear relationships between CT H and d C···M (R 2 = 0.901), ∆d av. MX (R 2 = 0.911) or D 0 (R 2 = 0.930) have been found for the analyzed (PH 3 ) 2 C· · · ZnX 2 dimers.  Figure 12. Relationships between d C···M and D 0 obtained for (PH 3 ) 2 C· · · MX 2 (X = Be, Mg, Zn) dimers.
The fact that it is imidazol-2-ylidene and (NH 3 ) 2 C molecules that willingly form different interactions with MX 2 than just C· · · M contact could have been expected from the electrostatic potential distribution map shown in Figure 5. From these maps, it is clear that both of these molecules have most acidic regions associated with the highly polar N-H bonds. It is therefore to be expected that these molecules will readily form a dihydrogen or hydrogen bond when the appropriate opportunity arises. Figure 13 shows that this is indeed the case. In cases (a) and (b), in addition to a magnesium bond to the carbon atom, there is also a dihydrogen bond. In the latter case it is very short, because the H· · · H distance is only 1.44 Å. It should be emphasized that it is not common for the dihydrogen bond between neutral molecules to be so short [191]. Of course, an acidic hydrogen atom from the N-H bond easily forms a hydrogen bond as well, as long as its acceptor is a strongly electronegative atom like fluorine, e.g., in MgF 2 (subfigures (c) and (d)). Again, this bond is clearly shorter (1.55 Å vs. 1.94 Å) when the N-H donor bond belongs to (NH 3 ) 2 C. In the case of imidazol-2-ylidene· · · ZnF 2 (e), the length of the N-H· · · F hydrogen bond is 1.92 Å, therefore similar to imidazol-2-ylidene· · · MgF 2 (c). However, in (NH 3 ) 2 C· · · ZnF 2 (f), the interaction between H and F is so strong that there is a highly advanced proton transfer to F, so that the distance H· · · F becomes much shorter (0.99 Å) than N· · · H (1.56 Å). Therefore, in this case, it is more logical to speak of a hydrogen bond of the F-H· · · N type. On the other hand, the dimers marked in Figure 13 as (g) and (h) are examples of systems with rather non-standard hydrogen bonds of the N-H· · · C type. In this pair, again, the interaction is much shorter for (NH 3 ) 2 C (1.87 Å) than for imidazol-2-ylidene (2.57 Å). The discussed examples are good illustrations of the coexistence of two formally completely different intermolecular interactions. Obviously, such an occurrence makes it much more difficult to extract the characteristic features for just one of them-in this case, the magnesium or zinc (spodium) bond.

QTAIM-and NCI-Based Characteristics
The characteristics of the studied dimers can be further investigated using the QTAIM [182][183][184] and NCI [192,193] theoretical methods. In particular, the former one is one of the most frequently used in studies of various intermolecular interactions. On the other hand, the latter of these methods is much less frequently used, and, to my knowledge, has not yet been utilized in the study of beryllium, magnesium, or zinc bonds.

QTAIM
As QTAIM has already been used previously for describing beryllium and magnesium bonds in some simple dimers [23][24][25][26][27], the main focus in this subsection is on the characteristics of the zinc bond and its possible differences from beryllium and magnesium bonds. For this purpose, QTAIM calculations were performed (ωB97X-D/6-311++G(2df,2pd)) for the following representative dimers: cyclopropenylidene· · · MX 2 (X = H, Br), imidazol-2-ylidene· · · MBr 2 , (PH 3 ) 2 C· · · MBr 2 , and (NH 3 ) 2 C· · · MBr 2 , where M = Be, Mg, Zn. These dimers were also chosen because they are examples of the dimers in which the previously described accompanying interactions either do not exist or do not have so significant influence on the C· · · M bond. Values of the most important quantities obtained by means of QTAIM are shown in Table 5. Table 5. Some fundamental QTAIM-based parameters (in au) characterizing carbene· · · MX 2 and CDP· · · MX 2 dimers: the electron density (ρ C···M ), its Laplacian (∇ 2 ρ C···M ) and the total electronic energy density (H C···M ) computed at the C· · · M bond critical point, the delocalization index of the C and M atomic basins (δ(C,M)). It is worth noting at the beginning that in terms of the obtained electron density (ρ C···M ) values or the total electronic energy density (H C···M ) calculated at the critical point of the C· · · M bond, the zinc bond does not differ much from the beryllium bond, whereas the corresponding values determined at the critical point of the magnesium bond are clearly different. For example, for the cyclopropenylidene· · · MH 2 dimer, the following values of ρ C···M and H C···M , respectively, for M = Zn, Be, and Mg have been obtained (in au): 0.074 ≈ 0.076 ≈ 0.033 and −0.023 ≈ −0.022 ≈ 0.003. Similarly, the corresponding pairs of triples of values for imidazol-2-ylidene· · · MBr 2 are: 0.093 ≈ 0.082 ≈ 0.046 and −0.032 ≈ −0.031 ≈ 0.000 and for (PH 3 ) 2 C· · · MBr 2 : 0.089 ≈ 0.085 ≈ 0.045 and −0.032 ≈ −0.036 ≈ −0.001. The similarity in terms of ρ C···M and H C···M of the zinc bond to the beryllium bond found here is an important result because, unlike the beryllium bonds [23][24][25][26][27][28][29][30], the former are studied only sporadically [91,93].

Dimer
For various types of interactions, the value of the electron density determined at the bond critical point of a given interaction (bond) is often treated as a measure of the strength of this interaction [182]. If so, then the zinc bonds should have a similar strength to the beryllium bonds (or even they should be slightly stronger than them), whereas the magnesium bonds should be much weaker. Comparison of the corresponding values of ρ C···M for the dimers with the same MBr 2 molecule suggests that cyclopropenylidene dimers should be the weakest, whereas (NH 3 Table 5. The corresponding relationships are shown in Figure 14 (left). They illustrate that, indeed, the relationships between ρ C···M and D 0 are very good (R 2 > 0.975), as long as the dimers with different M are treated separately.   Figure 14. The relationships between the electron density at the bond critical point of the C· · · M interaction (ρ C···M ) and either the dimer dissociation energy (left) or the C· · · M distance (right) determined for the dimers in Table 5.
It should also be interesting to check the quality of the linear relationship between the values of ρ C···M and d C···M . This relationship for the different M atoms is also shown in Figure 14 (right). For Mg and Zn, the coefficients of determination are very good (0.974 and 0.950, respectively), whereas the linear correlation is clearly worse for Be (0.888). This may result from much shorter C· · · M distances, and thus stronger intermolecular interactions. It is worth noting here similar slopes of the fitting lines for Be and Zn, which again supports the previously shown similarity of the beryllium and zinc bonds in the systems considered. On the other hand, the slope of the appropriate linear fit for Mg is much smaller, thus showing a much weaker relationship between the C· · · M distance and ρ C···M in the analyzed group of dimers.
The positive values (Table 5) of the Laplacian of the electron density at the bond critical point of C· · · M (∇ 2 ρ C···M ) show that this interaction is of closed-shell type [182]. However, all the complexes with Zn or Be taken into account in Table 5 feature a significantly negative value of the total electronic energy density at the bond critical point of C· · · M (i.e., H C···M ), which characterizes interactions with high degree of electron sharing, which in turn reflects a high degree of the C· · · M bond covalency [194]. On the other hand, the dimers with a magnesium atom feature H C···M values close to zero.
A very important QTAIM parameter often used to describe the A-B bond strength [195] is the so-called delocalization index δ (A,B), which defines the exchange of the electrons in the basins of atoms A and B [182][183][184]. Interestingly, particularly high values of δ(C,M) characterize the zinc bond, especially in (NH 3 ) 2 C· · · ZnBr 2 (0.732 au) and (PH 3 ) 2 C· · · ZnBr 2 (0.629 au), i.e., the systems in which CDP acts as the carbon atom donor. In a clear contrast, the δ(C,M) values for the dimers involving magnesium are similar to the δ(C,M) values for the dimers containing beryllium and are significantly lower than those for the dimers with zinc. Thus, surprisingly, the zinc bond to a carbene or CDP carbon atom should be much stronger than the corresponding magnesium or beryllium bond, of course, provided that the delocalization index is indeed a good measure of bond strength [195]. It is worth checking at this point whether there are strong linear relationships between the determined δ(C,M) values and other parameters describing the C· · · M bond strength, such as D 0 , d C···M and ρ C···M . These relationships are shown in Figure 15.   Figure 15. Relationships between the delocalization index of the C and M atomic basins (δ(C,M)) and the dimer dissociation energy (D 0 ), the electron density at the bond critical point of the C· · · M interaction (ρ C···M ), or the C· · · M bond distance (d C···M ) determined for the dimers in Table 5.
The quality of the obtained linear correlations clearly depends on both the correlated parameters and the type of the metal atom in the MX 2 molecule. In the case of the relationship between δ(C,M) and D 0 , the obtained linear correlations are reasonable for Mg and Zn (R 2 is 0.948 and 0.909, respectively), whereas the correlation for Be is rather weak (R 2 = 0.796). For the relationship between δ(C,M) and ρ C···M , the linear correlations are not great, especially for Zn (R 2 is ca. 0.9 for Be and Mg and 0.8 for Zn). In the case of the relationship δ(C,M) vs. d C···M the R 2 values for Be and Mg are pretty good (0.950 and 0.964, respectively), whereas for Zn, the linear correlation is clearly worse (R 2 0.865). It is worth noting that the obtained fitting lines for Zn are characterized by greater slopes, which of course results from the greater range of δ(C,M) values, from 0.493 au to 0.732 au (Table 5), thus ca. 0.24 au. In the case of Be and Mg, the range is only 0.08-0.10 au. This result indicates a greater number of electrons shared between C and Zn atomic basins than between C and either Mg or Be. Moreover, the amount is more dependent on the type of carbene or CDP.

NCI
Most of the QTAIM-based parameters are determined at critical points (e.g., of the C· · · M bond), and therefore these parameters are local, i.e., they provide information about the properties at a particular point in space. One way out of this limitation is the NCI method [192,193], which is based on the value of the reduced electron density gradient, s = 1/(2(3π 2 ) 1/3 )|∇ρ|/ρ 4/3 . Then, various interactions (especially those corresponding to low-density and low-gradient values) can be isolated by using appropriate cutoffs on the electron density values and its gradient. By means of the electron density gradient isosurfaces, individual interactions (especially non-covalent ones -hence the name of the method) show themselves as certain broad regions of real space rather than simply as a bond critical point between a pair of atoms [192]. In order to further investigate the difference between the zinc bond to the carbene or CDP carbon atom and its beryllium or magnesium counterpart, the electron density gradient isosurfaces were determined for the dimers of cyclopropenylidene, imidazol-2-ylidene, (PH 3 ) 2 C, and (NH 3 ) 2 C with MBr 2 , where M = Be, Mg and Zn ( Figure 16). The subfigure (a3) shows that the zinc bond in the cyclopropenylidene· · · ZnBr 2 dimer does not differ significantly from both Zn-Br bonds and should be stronger than the beryllium and magnesium bonds in its counterpart dimers with BeBr 2 (a1) and MgBr 2 (a2), respectively. It is worth noting that in the former of these cases, i.e., in the cyclopropenylidene· · · BeBr 2 dimer, two symmetrically located areas of weak interaction appear in the antibonding regions of the Be-Br bonds. This should lead to some elongation of both Be-Br bonds. The interaction picture in the case of imidazol-2-ylidene (b1-b3) is practically similar. In the case of the systems with MgBr 2 (b2) and BeBr 2 , (b1) one and two regions of very weak N-H· · · Br hydrogen bonds are visible, respectively, which are not followed by the presence of the respective bond paths. It has been shown that the presence or absence of a bond path generally has little to do with the interaction strength [30,196,197]. Furthermore, in the case of complexes with either (PH 3 ) 2 C or (NH 3 ) 2 C, representing the group of carbodiphosphoranes, the characteristics of changes in the areas of weak interactions caused by the change of the Zn atom to Mg or Be are similar. Specificially, in the case of the presence of the ZnBr 2 molecule (c3 and d3), the area for the C· · · Zn bond is similar to the area for both Zn-Br bonds, although it is distinguished by higher electron density values, especially in (NH 3 ) 2 C. In the (PH 3 ) 2 C· · · BeBr 2 (c1) and (NH 3 ) 2 C· · · MgBr 2 (d2) dimers, small areas of weaker interaction develop in the antibonding zones of the metal atom and in the (NH 3 ) 2 C· · · BeBr 2 (d1) dimer, these regions merge with the regions that characterize Be-Br bonds. These interactions, however, are clearly weaker than C· · · M, especially when the metal atom is Be. It can also be seen that in all the CDP-mediated dimers, there are two symmetrical P/N-H· · · Br hydrogen bonding regions, which are or are not (case c1) followed by bond paths. However, they should be much weaker than the C· · · Zn and C· · · Be bonds.
Summing up, it can be concluded that the analysis based on the NCI method shows that the zinc bond is the strongest, and although the beryllium bond should only be slightly weaker than it, the latter is related to the presence of additional areas of weaker interaction in the antibonding regions of the Be atom. The high strength of the C· · · Zn bond (competing even with the Zn-Br bond) is reflected in high values of ρ C···M and δ(C,M), but also in negative values of H C···M .

Theoretical Methods
Geometries of monomers and dimers were fully optimized on the ωB97X-D/6-311++G(2df,2p) level of theory, that is utilizing the ωB97X-D exchange-correlation functional [198] of Density Functional Theory (DFT) [199][200][201] and the 6-311++G(2df,2p) basis set [202][203][204][205][206], which includes both polarization and diffuse functions. By testing 200 different exchange-correlation functionals, the ωB97X-D functional has recently been shown [207] to be one of the best for general purposes. To increase the accuracy of the optimization procedure and numerical integration, cutoffs on forces and step size that are used to determine convergence were additionally tightened (0.000015 and 0.000010 for maximum force and its root mean square, respectively, and 0.000060 and 0.000040 for maximum displacement and its root mean square, respectively) and integration grid was increased to the (99, 590) one (UltraFine) having 99 radial shells and 590 angular points per shell. All the obtained systems were subjected to vibration analysis in order to check whether they correspond to the real minima on the potential energy hypersurface. There were no imaginary frequencies. Both geometry optimization and vibration analysis were performed by means of Gaussian 09 [208]. NBO-based [180,181] atomic charges were computed by means of NBO6.0 program [209] implemented in Gaussian 09. Calculations based on the QTAIM [182][183][184] and NCI [192] methods were made with the AIMAll program [210].

Conclusions
To date, the vast majority of theoretical studies on beryllium and magnesium bonds have used as Lewis bases small molecules, and the research on zinc (spodium) bonds is very rare. On the other hand, the research on carbenes and carbodiphosphoranes is mostly experimental. This article presents the results of theoretical research on the properties of beryllium, magnesium, and zinc bonds in a large group of dimers formed by the MX 2 molecule (where M = Be, Mg, Zn and X = H, F, Cl, Br, Me) and either carbene ((NH 2 ) 2 C, imidazol-2-ylidene, imidazolidin-2-ylidene, tetrahydropyrymid-2-ylidene, cyclopropenylidene) or carbodiphosphorane ((PH 3 ) 2 C, (NH 3 ) 2 C). Due to the rarity of theoretical studies of zinc bonds, the main focus in this article is placed on comparing them with both the beryllium bond and the magnesium bond.
The general characteristics of the presented dimers showed that the dissociation energies of the C(2)· · · M intermolecular interaction have wide range, from 10 to 53 kcal/mol, and this interaction is the strongest for the BeBr 2 and BeCl 2 Lewis acids. Although the C(0)· · · M bonds formed by (PH 3 ) 2 C are similar in strength to the C(2)· · · M bonds formed by carbenes, (NH 3 ) 2 C forms much stronger complexes, with a bond strength of up to 84 kcal/mol for the dimer with BeBr 2 . The interaction between MX 2 and either carbene or carbodiphosphorane leads to a significant bend of the MX 2 molecule, elongation of the MX bonds, and opening of the LCL angle (with a few exceptions). Importantly, it has been shown that the investigated systems are characterized by very high charge transfer effect from the carbene or carbodiphosphorane molecule to the MX 2 one. Even the weakest effect is more than twice as high as in the water dimer, while it is more than six times as strong in the (NH 3 ) 2 C· · · BeBr 2 dimer.
Theoretical studies based on the QTAIM and NCI methods have shown that the zinc bond is not very different from the beryllium bond; both should be of similar strength, while the magnesium bond should be weaker. Both are also characterized by a high degree of covalence. The determined values of the delocalization index show, however, that the zinc bond should be definitely stronger than the beryllium and magnesium bonds.
A large number of tested dimers as well as parameters characterizing both the interacting subsystems and the C· · · M bond itself allowed for the study of many linear relationships between the parameters. In general, they are good as long as systems with different M metal atoms are treated separately. The linear correlations for the zinc atom are usually slightly better than for the other atoms.
In addition to the dominant C· · · M interaction, some of the studied dimers also have various additional interactions, such as, e.g., the N-H· · · F, N-H· · · C and F-H· · · N hydrogen bonds, or N-H· · · H-Mg dihydrogen bond. In the latter case, it may be extremely short, such as 1.44 Å in (NH 3 ) 2 C· · · MgH 2 . These interactions, however, are much weaker than the beryllium, magnesium, and zinc bonds that are the main topic of the research.
A side result of the presented research is that the atomic charges obtained by the QTAIM method are highly unreliable. While more reliable than these, the NBO-based atomic charges also appear to be questionable. In contrast, the Hirshfeld atomic charges appear to be chemically sound.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.
Sample Availability: Optimized geometries of all considered systems are available from the author on request.

Abbreviations
The following abbreviations are used in this manuscript:

CDPs
Carbodiphosphoranes MESP Molecular electrostatic potential NBO Natural Bond Orbital QTAIM Quantum Theory of Atoms in Molecules NCI Noncovalent Interactions (method)