The many-body expansion for metals. I. The alkaline earth metals Be, Mg, and Ca

The purpose of the molecular many-body expansion (MBE) in quantum chemistry is twofold: in analysis, to understand the composition of interatomic and intermolecular interactions, and in practice, to facilitate calculations of large systems.^1–12 While computational chemistry has come far in both accuracy and efficiency, the scaling of most correlated methods is prohibitive for many practical applications. The MBE approach, on the other hand, relies on a series of smaller calculations of increasing n-body subsystems (i.e., dimers, trimer, tetramers) that aim to capture the most important properties and interactions of larger structures. Various MBE-based approaches have been extensively used to study many-water clusters,^13–17 water–ion interactions,^18,19 enzyme activity,²⁰ silicon nanowires,²¹ and much more. MBE-based methods are particularly attractive in situations when the many-body expansion converges rapidly, allowing for the MBE to be truncated to only the first few terms. This effectively reduces a large (sometimes impossible) calculation into multiple smaller calculations that can be easily parallelized. Recently, the capability of performant MBE-driven Molecular Dynamics (MBE-MD) simulations was also noted.²² Generally speaking, MBE methods are much more prevalent in non-covalently bound systems as non-covalent interactions can be straightforwardly decomposed into “bodies” and the expansion is rapidly convergent due to the long-range nature of their interactions.

The application of MBE in covalent structures, especially metallic systems, is fraught with challenges due to the complex electronic structure and short-range interactions exhibited by such systems. Both dynamic and non-dynamic correlation may be necessary to capture the multi-coordinated bonding of metals, in addition to considerations of different spins among the intermediate structures of open shell metals. For example, sodium clusters^23–26 (or any alkali metal) with an electronic configuration ns¹ will be in either a singlet or a doublet state at the various stages of the MBE (i.e., for the intermediate subsystems), and that is without considering the effect of additional states that might be relevant.^24,25 Transition metals have an even more complicated multi-state manifold^27–29 that could make identification of states challenging for intermediate terms on the many-body expansion.

The present article is the first in a series focusing on the application of the MBE in metallic systems. In this study, we focus on the alkaline earth metals (Be_n, Mg_n, Ca_n, n = 2, 3, 4, 5, 6). Alkaline earth metals have an even number of electrons with a closed-shell electronic configuration of ns² and primarily form singlet ground state clusters. Previous studies using density functional theory (DFT) methods have found high-spin (triplet and quintet) ground states for larger clusters of beryllium;^30,31 however, this is likely an artifact of DFT, as more accurate coupled-cluster calculations suggest a singlet ground state up to at least the beryllium hexamer.³² 

Although alkaline earth metals have a helium-like closed-shell electronic configuration, they exhibit markedly different bonding interactions.³³ Helium atoms interact primarily through weak van der Waals (vdW) forces^33,34 (noblesse oblige). Alkali earth compounds, on the other hand, resemble helium only in their dimer structures, which are also bound by weak vdW forces.³⁵ Larger alkali earth structures exhibit significantly stronger binding with large cohesive energies.^35–39 The large change from dimer to trimer has been primarily attributed to the stabilizing effect of non-additive many-body effect,^37,40 further explained as the result of the promotion of electrons from the s to the p orbitals in a multiconfigurational wavefunction, i.e., sp hybridization.^41,42

The literature on the binding of alkaline earth metals is quite extensive (see Refs. 41 and 42 and references within them); however, most studies seem to focus on the lowest energy (e.g., the tetrahedral) geometries of these structures. In order to practically employ MBE methods in the study of metallic systems (e.g., via MBE-driven molecular dynamics simulations), it is important to quantify and understand the behavior of the expansion in various geometrical conformations. In this paper, we examine several different geometries for the clusters of all three metals. Additionally, it is also worth investigating the binding of alkaline earth metals from a molecular orbital perspective, focusing particularly on the hybridized orbitals and their interactions. The outline of the paper is as follows: The computational methods used in this study are explained in detail in Sec. II. In Sec. III, the equilibrium geometries, binding and MBE energies, and the bonding analysis are presented and discussed. Conclusions are summarized in Sec. IV.

The geometries of the clusters examined in this study are shown in Fig. 1. All geometries were optimized at the Møller–Plesset second order perturbation theory (MP2)^43,44 using the aug-cc-pVDZ basis set.^45–47 To ascertain the effects of the different basis sets on the equilibrium cluster geometry, the tetramers were also optimized with the aug-cc-pVQZ basis set. The non-augmented version of these basis sets was used for the calcium clusters. The many-body energy terms were computed for all systems at the MP2 level of theory with the aug-cc-pVDZ and aug-cc-pVQZ^45–47 basis and the Coupled-cluster Singles and Doubles with perturbative Triples [CCSD(T)]⁴⁸ method using the aug-cc-pVDZ basis set. Only valence electrons were correlated with the MP2 and CCSD(T) methods. To test multi-reference effects on the MBE, calculations at the Multi Reference MP2 (MRMP2)^49,50 and the internally contracted Multi Reference Configuration Interaction (icMRCI)⁵¹ levels of theory were performed for a limited subset of the structures examined. The active space for the multi-reference calculations included all the valence s and p orbitals, totaling two electrons in four orbitals per atom (e.g., 8 electrons in 16 orbitals for the tetramers). MRMP2 calculations were also used to obtain potential energy surface (PES) scans of the dimers and the trimers of D_3h symmetry in order to analyze the evolution of the MBE along different bond lengths. The CCSD(T) calculations were performed using the NWCHEM quantum chemistry package,⁵² whereas all icMRCI calculations were performed with the MOLPRO 2015.1 package.^53,78 All other calculations were performed using the General Atomic and Molecular Electronic Structure System (GAMESS).^54–56 

There is unfortunately no silver bullet regarding the methods best suitable for MBE calculations of metals. MRCI might be the most accurate method; however, it is practically unfeasible for all but a few small clusters of light atoms; furthermore, it suffers from size-consistency and size-extensivity issues.^57,58 MRMP2 suffers from the same shortcomings, in addition to the possibility intruder states^51,59–62 (although this is rather unlikely in this particular case). CCSD(T), on the other hand, is computationally prohibitive for large systems. MP2 is perhaps more practical, though its accuracy remains to be validated. Both MP2 and CCSD(T) are single-reference methods, which may be inadequate in certain cases.

The MBE for the atomic clusters was calculated using the following formulas:^3,14,28,63

Higher order terms are similarly defined. For a system of N atoms, E_1B, E_2B, and E_3B, are the total one-body, two-body, and three-body energies, respectively, of the N-particle system. In Eq. (1), E_i is the total energy of atom i, E_ij in Eq. (2) is the total energy of the ij dimer, and so on. The dimer and trimer summations in Eqs. (2) and (3) run over all possible combinations of dimers and trimers, represented by the usual combination formula $CkN=N!k!N−k!$⁠.

The total binding energy, E_B, of the cluster with respect to the individual monomers (i.e., the separated atoms) is defined as

Finally, we have performed a quasi-atomic bonding orbital (QUAO) analysis^64–67,77 to understand the interactions of the MBE decomposition at the molecular orbital level. QUAOs are a type of localized orbitals with atomic orbital characteristics that can provide both qualitative and quantitative insight on the bonding patterns of molecules. The QUAOs are then used to construct the one-particle spin-averaged density,

In Eq. (5), the term p_Aa,Bb is the bond order (BO) density matrix computed in the basis of the QUAOs, denoted as $Aa$ (the notation means that orbital “a” is centered on atom “A”). The diagonal elements of the density matrix are the occupation numbers of the QUAO, while the interatomic off-diagonal terms represent the bond orders between QUAOs. Another useful metric is the kinetic bond order (KBO), defined as

The kinetic bond order correlates with the strength of the interaction between any two given QUAOs. It should be noted that the kinetic bond order does not necessarily correspond to the absolute binding energies between atoms—rather, it is a relative metric to compare the different bonding interactions between the QUAOs of the same molecule.

The QUAOs are computed at the full valence complete active space self-consistent field (CASSCF) level of theory, which includes all valence s and p orbitals and, in the case of calcium, also the d orbitals. Because of the computational cost and difficulty of full valence calculations, we focus only on a few key structures that help elucidate the binding of various MBE terms among these metals.

The geometrical parameters of the structures optimized for the beryllium, magnesium, and calcium clusters are listed in Table I. This collection of isomers is by no means exhaustive, but it represents a good sampling of one-, two-, and three-dimensional structures. Very often, only the lowest minima structures of alkaline earth metals have been examined in the literature. However, since our long-term goal is to employ the MBE in realistic applications (e.g., MBE-driven molecular dynamics), it is important to investigate its efficacy in structures other than global minima. Nevertheless, it was not possible to optimize all isomers for all three different metal clusters. For example, the zigzag geometry in beryllium always converged to either a linear or a pentagonal geometry, while the double triangle in magnesium converged to the hexagon geometry. This is reflected in the values presented in Table I.

Similarly, not all structures optimized correspond to local minima. For instance, for Be₄, only the tetrahedral and linear geometries are minima at the MP2/aug-cc-pVDZ level of theory, with the tetrahedral being the lowest energy structure. The rhombic and triangular geometries are single order saddle points, while the square is a fourth order saddle point. Interestingly, a minimum energy path⁶⁸ calculation of the square geometry leads to the tetrahedral geometry. Magnesium and calcium clusters show similar behavior, except that the square configuration is only a second order saddle point. The trigonal and tetragonal bipyramid configurations were the lowest energy minima for the pentamer and hexamer structures. Previous studies^30,31 have suggested that beryllium may have a high spin minimum, but in all our calculations, the singlet state was the lowest energy state at the MP2 or CCSD(T) levels of theory.

As expected, the increase in atomic radius down the group corresponds to a larger bond length (cf. Table I), although this increase is not uniform. For instance, the bond of the square structure increases by 1.44 Å from Be₄ to Mg₄ but only by 0.46 Å from Mg₄ to Ca₄. On average, the bond length increase from magnesium to calcium is half that of the increase from beryllium to magnesium. Conversely, the bond does not appear to change significantly as the cluster size increases from the tetramers to the pentamers and hexamers.

Table II compares the optimized geometries of the tetramers computed at different basis and levels of theory. As expected, the smaller double zeta quality basis set produces longer bond lengths, though in most cases, the differences are quite small. For instance, the bond distance of the tetrahedral structures (global minimum) changes by at most 0.1 Å between basis and levels of theory. Similarly, the difference between MP2 and numerically optimized CCSD(T) geometries is trivial in most cases. The few exceptions involve the square geometries, which differ by as much as 0.3 Å in bond lengths compared with the double zeta (DZ) optimized structures. This may be due to the strong multi-reference effects or unstable nature of the square geometry, which is located on a high order saddle point along the PES.

Tables III–V show the binding energies and MBE terms for all geometric conformations of beryllium, magnesium, and calcium clusters examined in this study. The convention that has been previously used in the literature^28,35,36,40 was to cast the energies normalized with respect to the two-body term, while preserving the original sign. Given that the three-body term is often the largest component of the MBE (vide infra) in such systems, we have instead chosen to normalize the MBE with respect to this largest three-body term. In this way, the MBE is more stable as an expansion as even small changes in the smaller two-body term due to electron correlation or the size of the basis set do not have a large effect on the ratios of the higher order terms. We, therefore, report the MBE as

A few properties of the MBE stand out immediately from Tables III–V, as well as Tables VI and VII, which compare the tetramer clusters at various levels of theory. For almost all cases and levels of theory, the non-additive three-body term is the largest in absolute energy and the most attractive term. Moreover, the change of energy between terms is often oscillatory: It decreases from two-body to three-body, increases from three-body to four-body, and so on. This is evident in Fig. 2, which shows a graphical representation of the MBE terms (in kcal/mol) for a few representative structures. Whether one considers an attractive or repulsive term in the MBE, beryllium always exhibits the largest magnitude in energy, followed by calcium and then magnesium. This order is also reflected in binding energies, in which calcium structures are often more strongly bound than magnesium (though weaker than beryllium). This energetic discrepancy in the trend between Be, Mg, and Ca is examined in Sec. III C. Nevertheless, it is uncertain if this trend continues in atoms below calcium, making beryllium a lone exception analogous to the first-row anomaly in the case of p-block elements,⁶⁹ or if a trend is absent altogether.

The question of whether the many-body expansion converges and at which rank is perhaps more complicated as different geometries exhibit radically different behavior in the MBE. The three-body term is unequivocally the largest term (in most geometries) and, therefore, must be included in any truncated application of the MBE. In most cases, the expansion appears to converge after the three-body term, though at substantially different rates (Tables III–V). For example, the repulsive four-body term in the tetrahedral structures of Be and Mg is smaller than the repulsive two-body term, but this is not so in the calcium tetrahedral geometry. The MBE of the trigonal bipyramidal structure converges rapidly for all atoms, while the pentagon is mostly non-convergent [at least at the CCSD(T) level of theory]. The oscillatory nature of the expansion also means that a truncated MBE can easily overestimate or underestimate the binding energy depending on where the truncation occurs. Therefore, the knowledge of the behavior of the MBE for a particular geometry is pivotal before the application of truncation. We plan to examine the application of truncated MBEs on alkaline earth metals in a future study.

Tables VI and VII show the comparison of the tetramer MBE between different computational methods [along with Tables III–V, which contain the CCSD(T)/aug-cc-pVDZ values]. There are several categories in which the data are compared: (1) MP2 vs CCSD(T), (2) aug-cc-pVDZ vs aug-cc-pVQZ basis, (3) results obtained with single- vs multi-reference methods, and (4) the effect of basis set superposition error (BSSE) corrections.

MP2 generally appears to overestimate the repulsive two-body and the attractive three-body terms (i.e., more positive two-body energies and more negative three-body energies) compared to CCSD(T). This is consistent with previous observations in the literature.³² For example, the MP2 three-body term of the Be₄ square geometry is 5.8 kcal/mol larger than the corresponding CCSD(T) energy, while the MP2 three-body term for the Be₄ tetrahedral geometry is 53.7 kcal/mol larger, and so on. Similar patterns are observed for all magnesium and calcium structures. Because the depictions of the MBE are normalized with respect to the three-body terms, the normalized values remain relatively stable. In general, a larger two-body energy also corresponds to a larger three-body energy, a fact that maintains a proportional change among the MBE terms. This is graphically demonstrated in Fig. 3, which shows the correlation between the energy differences between the two-, three-, and four-body terms with respect to the (often larger) three-body energies. It is evident from the fitted lines that the relationship between the MBE terms does not change much between the MP2 and the CCSD(T) calculated values. Both methods are fitted to similar equations and accuracy (R² value). In fact, even the available icMRCI data (red and green points in the right graph of Fig. 3) follow the fitted lines quite well. This suggests that regardless of the apparent differences between computational methods, the qualitative behavior of the MBE (such as the ratio of the terms with respect to the three-body term) remains the same.

On the other hand, the choice of basis set (Table VI) appears to have a larger and more complicated effect on the MBE. For most geometries optimized with the aug-cc-pVDZ basis set, the MBEs obtained with the larger quadruple-zeta (aug-cc-pVQZ) basis set show a decrease in the two-body energy while maintaining the three- and four-body energies relatively constant (for a comparison of unscaled MBE energies, see Table S2 and Figure S4). For example, in the rhombic beryllium tetramer, the repulsive two-body energy is nearly halved in the larger basis set, while the three-body and four-body energies only change by a few kcal/mol (cf. Table VI). In the tetrahedral geometries, for example, the three-body energies in the QZ basis decrease only by 3.2 kcal/mol, 0.5 kcal/mol, and 1.5 kcal/mol for Be, Mg, and Ca respectively. Computing the MBE with the aug-cc-pVQZ basis set for geometries also optimized with the aug-cc-pVQZ basis set has a more significant effect, resulting in a slight decrease in both the two- and three-body terms. Although in most geometries, these changes are small enough to be qualitatively similar, there are some substantial outliers (cf. Table VI). The three-body energy of the tetrahedral geometry, for example, decreases by ∼17 kcal/mol with the aug-cc-pVQZ basis set at the aug-cc-pVQZ geometry, considerably more than the changes in the two-body or three-body energy (or the changes in any other geometry), despite the very small change in bond length (see Table II). This is likely because the tetrahedral structure has the shortest interatomic distances of all the geometries and, therefore, even a small change in geometry can affect the bonding interactions of the hybridized orbitals (discussed in Sec. III C).

It is clear that there is a basis set dependence of the MBE. This effect is larger for the two-body energy terms because they are generally smaller in energy than the three-body and four-body terms (Fig. 2). Thus, even a small change in energy can be a significant portion of the total two-body energy. In geometries that are weakly bound (or outright repulsive), such as the linear and square geometries of magnesium and calcium, the two-body term can be the largest portion of the MBE (see Table VI); therefore, it can be significantly affected by the choice of basis set. In fact, the geometrical differences of these structures can be explained by the dominance of the attractive two-body term, which causes the reduction of the bond lengths. Table II shows that geometries optimized with the CCSD(T)/aug-cc-pVQZ level of theory can result in as much as a 0.3 Å reduction in the bond length. Considering the MBE of some of these structures at the CCSD(T)/aug-cc-pVQZ level of theory,

Despite the significant changes in the bond length, the CCSD(T)/aug-cc-pVQZ MBE follows the same trend as the MBE computed with the smaller basis set, albeit with a more attractive two-body term. The three-body terms also become more attractive, no doubt a result of the shorter bond lengths. Nevertheless, we emphasize that the overall trend of the MBE remains qualitatively the same in most cases (see the supplementary material for more examples).

Even at the atomic level, alkaline earth atoms exhibit some multi-reference character due to the symmetry allowed s–p double excitations, and this effect increases for the larger clusters, although this increase is not uniform among the different metals or the studied geometries. For instance, the highest occupied natural orbital occupation numbers of the beryllium, magnesium, and calcium atoms are 1.83, 1.85, and 1.83, respectively. In the D_3h trimer geometry, this occupation becomes 1.77, 1.83, 1.79 for Be, Mg, and Ca respectively (note that these are doubly degenerate orbitals of E symmetry). The linear tetramer geometry becomes more strongly multi-reference, with highest orbital occupation numbers of 1.11, 1.84, and 1.22, respectively. It is unclear what causes the smaller multi-reference character of magnesium compared to the other alkaline earth metals, as a visual inspection of the orbitals shows no appreciable differences between the three metals. It is likely that the larger core–core repulsion in magnesium dominates over the bonding interactions of the 2s valence electrons. Nevertheless, despite the presence of multi-reference character for many of the probed structures, the difference in the MBE composition between single- and multi-reference methods appears to be minimal (see Table VII). The CCSD(T), MRMP2, and icMRCI results follow the same qualitative trend. While the exact energies differ [icMRCI appears to be closer to CCSD(T) than MRMP2], the changes in energies are relatively similar. This is also shown in Fig. 3, where the available icMRCI calculations are superimposed on the CCSD(T) line. In many cases, the MRMP2 method appears to produce more attractive three-body terms (and similarly, more repulsive two-body terms), whereas the icMRCI results appear remarkably close to those of CCSD(T). For example, the beryllium tetrahedral three-body energy is −112.4 kcal/mol, −123.4 kcal/mol, and −113.4 kcal/mol for CCSD(T), MRMP2, and icMRCI, respectively. Similar trends are observed in the other geometries, except for the extended triangular geometry in which the MRMP2 produces a less attractive three-body energy (by ∼7 kcal/mol). Nevertheless, the differences are still qualitatively minimal. This is even more the case in magnesium in which there’s a smaller difference between the three methods (Table VII).

Finally, Table VIII shows the BSSE corrected energies for the MP2/aug-cc-pVDZ and MP2/aug-cc-pVQZ levels of theory. Counterpoise corrected MBE energies differ by as much as 1 kcal/mol from their uncorrected counterparts (naturally more for the aug-cc-pVDZ than for the aug-cc-pVQZ basis set), though the BSSE correction is not significant enough to make a qualitative difference in the normalized many-body terms.

Beryllium is unique among the alkaline earth metals because its overall binding energy is known to be stronger than that of magnesium and calcium,^33,39 a result that was also confirmed by our calculations. For instance, the binding energies of tetrahedral (Be₄), trigonal bipyramidal (Be₅), and tetragonal bipyramidal (Be₆) are −70.14 kcal/mol, −99.31 kcal/mol, and −104.2 kcal/mol, respectively. These are two to three times stronger than their magnesium or calcium counterparts (see Tables III–V). Figure 4 shows the potential energy surfaces of the dimers and the D_3h trimers (symmetric stretch) at the MRMP2/aug-cc-pVQZ level of theory (except for Ca, which used the cc-pVDZ basis) for all three metals examined here. The beryllium dimer is weakly bound with a binding energy of −974 cm⁻¹ (∼2.8 kcal/mol, including zero-point energy) at an equilibrium distance of 2.42 Å. These values compare reasonably well with previous experimental measurements (energy of 750–930 cm⁻¹ and distance of 2.45–2.54 Å, depending on the experiment).⁷⁶ The binding energy increases nearly tenfold for the trimer and even more so for the tetramer (not shown in the graph). This trend holds true for magnesium and calcium as well, though the trimers of those systems are not as strongly bound. For instance, the Mg trimer has a binding energy of 4.8 kcal/mol (8-fold increase from the dimer energy), while the Ca trimer has a binding energy of 6.15 kcal/mol (20-fold increase from the dimer energy). Beryllium clearly shows a much stronger three-body interaction compared to the metals below it in the Periodic Table, followed by calcium and then magnesium.

Similarly, Figs. 4(a)–4(c) document the repulsiveness of the two-body term at the corresponding trimer equilibrium distances. As the dimer minimum is located at larger metal–metal distances than the trimer (or tetramer, not shown in the graph), at the trimer equilibrium separation the two-body energy falls along the repulsive wall of the dimer PES. In the case of calcium, the two-body energy is almost 0 kcal/mol at the trimer equilibrium separation. This effect has been previously reported for the sodium trimer.^24,25

To compare the differences between the PES of the clusters, we examine the potential energy surfaces in reduced forms (expressed as a ratio with respect to the distances and energies at the equilibrium geometry of each respective structure)^70,71 shown in the second row of Figs. 4(d)–4(f). In these reduced coordinate forms, there is little difference between the dimer and total two-body surfaces of beryllium, magnesium, and calcium clusters. This is evident in the reduced coordinate surfaces [Fig. 4(d)], which shows surfaces virtually overlapping at equilibrium and long-range distances. The dimers of all three systems are bound equally weakly, a factor that explains their similarity in the reduced plots. In the short-range regions of the dimers, there is a clear progression from Be to Mg to Ca in which the energy becomes more repulsive down the Periodic Table. This behavior is typical of the increase in core–core repulsion and atomic radius with atomic number, analogous to the differences of ionic radii that have been previously reported.⁷¹ Note that for the D_3h trimer, the two-body energy is exactly three times that of the dimer; therefore, the reduced two-body curve would overlap that of the dimer.

On the other hand, the trimer reduced PESs are markedly different from each other without exhibiting much overlap between the three metals. The reduced three-body terms of magnesium and calcium are fairly closer together, while the reduced three-body energy of beryllium is noticeably smaller (i.e., more attractive). This is reflective of the stronger three-body term of beryllium. The short-range region of the curves follows the same pattern, with magnesium showing a higher repulsion than calcium and beryllium. We speculate that the presence of virtual d-orbitals in the valence space of the calcium atom enables it to form stronger attractive interactions than magnesium, counteracting the strong core–core repulsion (discussed below). Unlike the dimer curves, there is no a priori reason that the reduced trimer or three-body curves will overlap in the same manner as they do reflect different behavior among the metals.

A better understanding of these interactions from a molecular orbital perspective can be achieved via a quasi-atomic orbital (QUAO) analysis. Figures 5–11 display the most important QUAOs for the dimer, trimer, and tetramer of the tetrahedral geometry for Be Mg, and Ca metals, i.e., all interatomic distances are those of the respective tetrahedral arrangements. In addition to the density plots of each respective QUAO (shapes), several useful quantitative values emerge from the density matrix in the QUAO basis: the orbital occupation number (Occ), the bond order (BO), the kinetic bond order (KBO), and the atomic-orbital hybridization decomposition (percentage of s, p, or d mixing in the orbitals). The naming convention used in the figures below is as follows: QUAOs centered on atom X and directed toward atom Y are labeled as XY-c, where “c” represents the character of the orbital, such as s, p, d, sp (hybridized s and p), and lp (lone pairs). For example, an sp orbital centered on the one beryllium atom (Be₁) and directed toward another beryllium atom (Be₂) is denoted as Be₁Be₂-sp. Conversely, if an orbital is centered on an atom X but not directed toward another, it is simply labeled as X-c. Note that the labels used are merely a descriptive convenience—that is, a QUAO labeled “lp” (lone pair) is not necessarily a doubly occupied nonbonding orbital.

The quasi-atomic orbitals for beryllium reveal an increasing degree of hybridization with increasing size of the cluster (Fig. 5). The beryllium dimer exhibits an approximate 50/50 mixture of s and p orbitals (sp hybridization). This increases to approximately an averaged sp² hybridization for the trimer and an averaged sp³ hybridization for the tetramer. The increase in hybridization certainly allows for more interactions in the larger structures (and stronger binding energies), but it may not be the only cause for the larger binding energies. It has been suggested in the literature that the sp hybridization is responsible for the binding of beryllium clusters.^41,42 The occupation numbers of the hybridized orbitals shown in Fig. 5 also support this finding as the occupation number seems to correlate significantly with the degree of hybridization. The pure p orbitals of the dimer and trimer are virtually unoccupied, while the sp-hybridized ones have large occupation numbers.

Examining the interactions between the beryllium QUAOs in Fig. 6 shows that the dimer is primarily bound by the interaction of two sp QUAOs and a much weaker lp–lp interaction. It is tempting to consider the sp–sp interaction as a σ bond; but, as shown in Fig. 4, this is a much weaker interaction than a conventional sigma bond of, say, a carbon atom. Similarly, there is no interaction whatsoever between the p orbitals perpendicular to the diatomic bond.

Unlike the dimer, the beryllium trimer and tetramer benefit from the availability of more interactions between the QUAOs. Figure 6 shows not only the direct sp–sp interaction between the Ca₁Ca₂-sp and Ca₂Ca₁-sp orbitals (one of which also exists in the dimer) but also a significant interaction between orbitals not directed toward their respective centers, viz., the interaction between Ca₁Ca₃-sp and Ca₂-lp orbitals. The second interaction is inevitably absent from the dimer and important to the stability of the trimer. Though energetically weaker than the sp–sp interaction by the KBO metric (−90.3 kcal/mol vs 2 × −15.7 kcal/mol per atom pair), it contains most of the bond order for the molecule (0.69 vs2 × 0.47 per atom pair). The same holds true for the tetramer, in which the sp–lp interactions contain most of the bond order. Therefore, it is very likely that the presence of these indirect sp–lp interactions, more so than the sp–sp interactions, contributes to the stability in the trimer, which, as Fig. 4 suggests, is much lower in energy than that of the dimer. The QUAO interactions in the trimer also help demonstrate the non-additivity of three- and higher-body terms. The sp–lp interaction between two beryllium atoms can only exist in the presence of a third atom.

On the other hand, the magnesium QUAOs tell a different story. As Fig. 7 clearly shows, the hybridization of magnesium is considerably smaller than that of beryllium. While the dimer of the two elements is more or less the same (sp-hybridized), the trimer and tetramer structures of magnesium show virtually no hybridization whatsoever. For example, in the trimer, the highest occupied QUAO (1.37 occupation number) has 96% s character and only 4% p character. Moreover, the virtually pure p orbitals have an occupation of 0.586, which is non-negligible but considerably lower than the hybridized orbitals for beryllium. Instead, most of the electron occupation is in the s orbital. The tetramer exhibits similar properties as the trimer, but with the p orbitals having a larger share of the occupation number.

The interactions between the QUAOs in the magnesium structures (Fig. 8) are also different from those in beryllium. Since there is not any appreciable sp hybridization in the trimer and tetramer structures, there are not any sp–sp interactions or sp–lp interactions equivalent to those of Be. Additionally, the s–s interaction is negligible. Instead, the predominant interactions are of s–p and adjacent p–p types, which also demonstrate the non-additive nature of the trimer interactions (though the energy difference between the Mg dimer and trimer are considerably smaller than those of Be).

The difference between magnesium and beryllium encountered here can be explained, at least partially, by the competition between the binding interaction of the valence electrons and the core–core repulsion. Although both Be and Mg (as well as Ca, discussed below) have two valence electrons, beryllium has only two core electrons, while magnesium has ten core electrons, allowing for a stronger repulsion between the atoms. This is evident from the increase in bond length of Mg compared to Be. The fact that the atomic radii and the average metal–metal bond length increase down the Periodic Table is well known,⁷² and in the case of magnesium, the electronic repulsion overcomes the valence–valence interactions. In beryllium, the repulsion of the two core electrons is not enough to counterbalance the valence–valence interactions.

The distance dependence of the interatomic interactions becomes evident once the QUAOs are computed at various points in the PES, as shown in Figs. 12 and 13. In the case of beryllium, even a slight increase of bond length leads to a loss of hybridization. Conversely, reducing the distance of Mg₃ introduces hybridized orbitals. In other words, at comparable distances both Be and Mg show considerable sp hybridization, which, in principle, leads to similar valence–valence interactions. However, the larger core and electronic repulsion of Mg pushes the atoms further apart, thus reducing the bonding interactions in comparison to beryllium. Note that optimizing the geometries with larger basis sets or different computational methods can result in slightly shorter bond lengths for the clusters examined here (Table II). In most cases, however, these changes are not enough to significantly affect the binding between atoms (such as the changes in the Mg or Ca structures). It is possible that the interactions between QUAOs can be affected appreciably with large changes of the atomic distances as a result of using larger basis sets (see Table II).

Computing the QUAOs for the square geometry (Fig. 9), one arrives at a similar conclusion. The bond distance of the square geometry is longer than that of the tetrahedral geometry (in the case of Be, it is longer by 0.7 Å at the aug-cc-pVDZ basis set, and 0.55 Å at the aug-cc-pVQZ basis). The QUAOs between the two systems are virtually the same (no sp hybridization), a fact also mirrored by comparable binding energies [−4.30 kcal/mol and −3.04 kcal/mol for Be and Mg, respectively, at the CCSD(T)/aug-cc-pVDZ level of theory]. Indeed, this behavior is quite revealing, as the distance dependence of sp hybridization appears to correlate strongly with the increased attractiveness of the three-body energy of the trimer (cf. Fig. 4). As the sp hybridization increases at shorter distances (Figs. 11 and 12), so does the three-body energy. This is a strong indicator of the importance of hybridization on the stability of these clusters.

Initially, one would expect calcium to follow the same trend as magnesium. Indeed, calcium has a larger atomic radius and stronger core–core repulsion, leading to larger equilibrium distances (see Table I and Fig. 4). However, the binding energy of calcium calculated here and in the literature³⁶ is larger than that of magnesium for all geometries and levels of theory. The main difference between calcium and beryllium or magnesium is the presence of low-lying d orbitals, which can actively participate in the bonding interactions (Ca has a 4s² configuration). This poses a significant computational challenge for the purpose of computing the QUAOs. Including d orbitals in the CASSCF calculation would require an active space of nine orbitals per atom (1s + 3p + 5d), which would be computationally feasible only for the small clusters. To sidestep this issue, we used the Occupation Restricted Multiple Active Spaces SCF method (ORMAS-SCF)^73,74 in order to partition the active space into two subspaces: one containing the s and p orbitals and another containing the virtual d orbitals. All possible configurations are generated in the s–p subspace (i.e., a full CI), but only single and double excitations are allowed in the d subspace. Double excitations into the d space should be sufficient to capture the contributions of the d orbitals into the interatomic interactions.

The calcium QUAOs shown in Fig. 10 indicate considerable mixing of the d orbitals with the s and p orbitals, especially in the trimer. The orbitals in the calcium dimer look no different than those of beryllium, except for a small amount of d mixing in the Ca₁Ca₂-sp orbital. The highest occupied orbital of the trimer structure, Ca₁-s, has as much as ∼10% d-character (i.e., more than p-character), while the Ca₁Ca₂-πd orbital has 21% d-character and an occupation of 0.30 electrons. The presence of d orbitals in the space of QUAOs likely “expands” their range and, unlike Be or Mg, enables interactions to occur at larger distances.

Figure 11 shows the most important interactions (measured by BO and KBO) between the QUAOs of calcium. Again, there is little difference between the calcium, beryllium, or magnesium dimers. However, the difference in the trimers is considerable. Unlike magnesium, there is a small s–s type interaction between atoms, likely due to the presence of d orbitals in the Ca₁-s QUAOs (the Mg₁-s QUAO had 96% s character). The most dominant interactions are s–p (BO = 0.59), s–πd (BO = 0.53), and p–πd (BO = 0.27). These interactions are also non-additive in nature (they are composed of orbitals present only in the trimer) and serve to stabilize the three-body energies of the trimer.

Unfortunately, the KBOs used here are unable to quantitatively distinguish between the QUAO interactions of different atoms (e.g., Mg vs Ca), but they do so for interactions only within the same cluster. However, the presence of d orbitals in the valence hybridization of calcium allows for more interactions than magnesium. Thus, the d orbitals are unequivocally involved in the bonding of the calcium trimer (unlike beryllium or magnesium), likely explaining the larger binding energy of calcium clusters.

It has been well known in the literature that beryllium structures have a stronger binding energy than magnesium and calcium, and the QUAO analysis presented here links this property directly to the strong hybridization in conjunction with the smaller core–core repulsion of the beryllium atom. However, this raises the question about the potential presence of atoms in their excited states in the cluster, a fact that will affect the MBE. So far, all energies needed for the MBE have been calculated strictly with respect to the ground states. While these clusters do dissociate to the atomic ¹S ground states (see PES graph in the supplementary material, Fig. 1), work by Kalemos⁷⁵ makes the case for the presence of ³P states at the equilibrium geometries of the beryllium dimer and trimer, based on the results of the atomic populations. Indeed, our QUAO analysis shows the significant presence of sp hybridization that increases with increasing cluster size and decreasing bond lengths (see Figs. 5, 12 and S3). This could be the result of the increasing multi-reference character of beryllium structures rather than an excited atomic state, though the precise determination of atomic states merits a separate study of its own and will not be examined here. Nevertheless, it is important to mention this factor, and for completeness, we present a few alternative MBE calculations of the beryllium trimer in Fig. 14.

It becomes clear from Fig. 14 that by considering intermediates in the MBE where the atoms are in their excited states, the determination of the MBE becomes increasingly complicated. In addition to the many possible combinations of states, one must also ensure that the correct selection rules and mixing of terms are observed. For instance, a $Σg+1$ can result from ¹S+¹S or ³P+³P states, but not from a ¹S + ³P combination. Additionally, both ¹S → ³P and $Σg+1→3Σu+$ transitions are spin forbidden (though dipole allowed), but they may be allowed under the influence of a perturbation that breaks the symmetry (i.e., the presence of another atom).

In accordance with Eqs. (1)–(3), the diagrams in Fig. 14 indicate the states of the monomers and dimers used to compute each two-body and three-body terms. While the dimer with respect to the ground state is repulsive, when considered with respect to the ³P asymptotic states of the atoms, it is attractive by −121.4 kcal/mol. However, considering a two-body term composed of ³P states results in a repulsive three-body energy (as the equation for the three-body contains the energies of monomers and dimers). In other words, a repulsive two-body term is required to obtain an attractive three-body term, with the magnitude of the energy depending on the state of the monomers—i.e., more ³P monomers correlate with a lower three-body energy.

Conversely, one may consider the first excited dimer of state $Σu+3$⁠, which lies 19.5 kcal/mol above the dimer ground state. This state dissociates to the ¹S + ³P atomic states, though such a state can also result from a ³P + ³P atomic combination (this is not considered in Fig. 14). This dimer state results in two attractive and one repulsive three-body term, with the latter deriving from ground state monomers.

For convenience, the various schemes of calculating the two- and three-body terms of the D_3h Be trimer are summarized in Table IX. Note that the binding energies for each of these schemes are the same (minus rounding errors) as they all start and end at the same states. The beryllium trimer is in the ¹A_g state at equilibrium geometry and dissociates into the three ¹S states; therefore, that is the correct binding energy for the system. It is obvious how radically the terms of the MBE can change when considering excited intermediate states, and how this could affect the convergence of the MBE.

In this study, we examined the interactions and the MBE of three alkaline earth metals, beryllium, magnesium, and calcium, for several geometry conformations and levels of theory. Our results confirm earlier reports that the non-additive three-body term is the largest and most important attractive term, followed by oscillatory and decreasing in magnitude higher order terms. However, the behavior of the MBE differed significantly among the geometries tested, showing different rates of convergence or non-convergence at all in some cases. Our results indicate that the qualitative trends of the MBE, such as the energy differences from one term to another, do not change significantly between levels of theory, with the biggest effect arising from the basis set, which in some geometries resulted in shorter bond lengths. Though these differences were not large enough to change the overall trend of the MBE, significant reductions in the bond lengths can lead to more attractive two-body or three-body terms. The inclusion of dynamic correlation was also important in being able to locate a minimum on the dimer PES in all three metals. Similarly, single reference methods appear adequate at describing the MBE, even though the multi-reference character of these systems was found to be non-negligible.

A detailed analysis of the dimers, trimers, and tetramers using quasi-atomic molecular orbitals helps provide an understanding of the many-body interactions in alkaline earth metals. As previously speculated, the QUAO analysis reveals that sp hybridization is crucial at stabilizing the three-body term of beryllium. This is absent in magnesium at equilibrium distances due to the large core–core repulsion forces that extend the average bond length compared to beryllium. This is also true for calcium; however, the presence of d-orbitals (in the form of hybridized s–p–d orbitals) facilitate stronger interactions in the Ca clusters.

Finally, the possibility of excited state intermediates in the many-body expansion of beryllium is also briefly considered, stemming from the sp hybridization of the beryllium structures. Different monomer and dimer states can lead to radically different three-body energies, and this certainly complicates MBE as it is far from trivial to determine which states may be relevant in a given term of the expansion. However, a more detailed examination of this issue is beyond the scope of this article and will be examined in a future study.

All structures optimized in this study are included in the supplementary material in xyz format. Also see the supplementary material for additional potential energy surfaces of the Be and Mg dimer and for population analyses.

J.M. was supported by the US Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences at Pacific Northwest National Laboratory (PNNL), a multi-program national laboratory operated for DOE by Battelle. S.S.X. was supported by the Center for Scalable Predictive methods for Excitations and Correlated phenomena (SPEC), which is funded by the US Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences and Biosciences Division as part of the Computational Chemical Sciences (CCS) program at Pacific Northwest National Laboratory. This research also used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the US Department of Energy under Contract No. DE-AC02-05CH11231.

The authors have no conflicts to disclose.

Joani Mato: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Demeter Tzeli: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Sotiris S. Xantheas: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.