Reply to the ‘Comment on “The oxidation state in low-valent beryllium and magnesium compounds”’ by S. Pan and G. Frenking, Chem. Sci., 2022, 13, DOI: 10.1039/D2SC04231B

A recent article by Pan and Frenking challenges our assignment of the oxidation state of low valent group 2 compounds. With this reply, we show that our assignment of Be(+2) and Mg(+2) oxidation states in Be(cAACDip)2 and Mg(cAACDip)2 is fully consistent with our data. Some of the arguments exposed by Pan and Frenking were based on visual inspection of our figures, rather than a thorough numerical analysis. We discuss with numerical proof that some of the statements made by the authors concerning our reported data are erroneous. In addition, we provide further evidence that the criterion of the lowest orbital interaction energy in the energy decomposition analysis (EDA) method is unsuitable as a general tool to assess the valence state of the fragments. Other indicators based on natural orbitals for chemical valence (NOCV) deliver a more reliable bonding picture. We also emphasize the importance of using stable wavefunctions for any kind of analysis, including EDA.

Pan and Frenking 1 have challenged our recent assignment of the formal oxidation state of the beryllium and magnesium family of di-coordinated compounds ML 2 (M = Be, Mg) with ligands of cyclic (alkyl) (amino) carbene (cAAC) and N-heterocyclic carbene (NHC) types. 2 In a landmark study, the group of Braunschweig isolated and characterized the Be(cAAC Dip ) 2 complex. 3 The authors described the system as the rst example of a stable beryllium compound in the zero oxidation state Be(0), where the bonding situation consists of two donor-acceptor interactions between cAAC ligands and the central atom cAAC(0) $ Be(0) % cAAC(0). The preparation of the heavier analogue Mg(cAAC Dip ) 2 has been attempted by Turner, but under the reaction conditions a ligand activation is observed, giving evidence for a highly reactive transient species. 4 In our original work, 2 we provided an alternative picture where the strong p-acidity of the ligand causes its own reduction giving rise to the singlet diradicaloid L(−1) / E(+2) ) L(−1) situation. 2 The non-innocent behavior of cAAC ligands in main group and transition metal complexes has been discussed in a number of computational and experimental studies. [5][6][7][8][9] One of the most relevant ndings of our work 2 is that the closed-shell single determinant wavefunction is either not stable (B3LYP) or higher in energy than a broken-symmetry KS-DFT solution. This indication of signicant static correlation prompted us to rely on multireference wavefunctions (CASSCF in particular) to study the electronic structure of these species. Besides the usual electronic structure indicators, we applied the effective oxidation state (EOS) analysis to assign the oxidation state of the central atom and that of the ligands. This approach, introduced in 2015 by some of us, 10 has already been tested on a wide variety of challenging systems, exposing some limitations of algorithms preferred by IUPAC. 11 In contrast to energy decomposition analysis (EDA), [12][13][14] EOS was specically devised solely for the purpose of assigning oxidation states. It should be noted that other approaches specically designed for OS assignment such as Head-Gordon's Localized Orbital Bonding Analysis (LOBA) 15 or the recently introduced Oxidation State Localized Orbitals (OSLO) 16 can only be applied for singledeterminant wavefunctions. EOS is the only scheme of that kind that can be applied on equal footing for single-determinant and multicongurational wavefunctions, and hence it was the natural choice for our study.
Remarkably, only a small part of the criticism of Pan and Frenking actually refers to our preferred methodological a Departament de Química, Institut de Química Computacional i Catàlisi, Universitat de Girona, c/M. Aurelia Capmany 69, 17003 Girona, Spain. E-mail: pedro. salvador@udg.edu framework, with arguments that are rather bizarre or simply faulty. For instance, they state that the EOS scheme "[.] is based on atomic fragments that are dened in different ways and can lead to different results". 1 While this is a pretty accurate description of the EDA approach, it does not apply to EOS analysis. In EOS, one needs to dene the composing fragments, in this case the central atom and each ligand, but not their reference state. As such, the EOS approach is free of bias (vide infra). 6 Then, the occupation of the fragment's domain natural orbitals (i.e. Mayer's effective fragment orbitals, 17 henceforth EFOs) is used to assign electron pairs (or individual electrons in open-shell systems) to the fragments. Rather than rounding the occupation numbers to the closest integer, the EFOs of all fragments are sorted in descending order. The most plausible OS assignment (for all fragments) is the one in which the electron pairs are assigned to the highest occupied EFOs. A reliability index R (%) quanties the extent to which the discrete OS values match the actual electronic structure. It is simply derived from the occupation numbers of the frontier EFOs in the EOS procedure. 10 In fact, it is not the way the fragments are dened, but rather the identication of the atom in the molecule (AIM) that inuences the numerical results to some extent. In Tables S9 and S10 of our original work, we already provided numerical evidence of the effect of using one or another AIM denition to obtain the EFOs. 2 Pan and Frenking pinpoint the case of Be(cAAC Me ) 2 , for which using classical Löwdin analysis results in a Be(0) OS (with a very low R (%) value, a close-call situation, in contrast with a clear Be(+2) OS when using more reliable AIM approaches like the Quantum Theory of Atoms in Molecules (QTAIM) 18 or the Natural Atomic Orbital (NAO) basis. 19 Pan and Frenking failed to notice that Löwdin's partial atomic charge on Be is −0.47 e, in sharp contrast with the other AIM approaches tested (well over +1.0 e). They also missed the case of Mg(cAAC Me ) 2 . Here, the Löwdin partial charge on Mg is as low as −0.75 e, while it is ca. +0.4 e when using NAO or QTAIM. Remarkably, in all three cases the EOS analysis still yields the same undisputable OS assignment of Mg(0). A less biased perception of Tables S9 and S10 indicates that (i) EOS analysis is signicantly more robust than partial atomic charges and (ii) partial atomic charges neither correspond to nor correlate with  oxidation states, as repeatedly observed, and in line with IUPAC's denition of OS and the winner-takes-all principle. 20 At this point it should become clear that we did not "omit the possibility that the metal atoms has the oxidation state +1": such a faddy result (which necessarily implies a formal OS of − 1 2 for each of the ligands) simply cannot stem from the wavefunction analysis. Articially forcing the EOS scheme to produce a particular set of OS on the fragments, by denition, results in a lower R (%) value (in fact, even below the worst-case scenario of R (%) = 50).
Pan and Frenking also cast doubts on the fact that for Mg(cAAC Me ) 2 we obtain Mg(0) while for Mg(cAAC Dip ) 2 it is Mg(+2). Accordingly, "a drastic alteration in the nature of the Mg-cAAC R bond from Mg(0) to Mg(2+) by changing the rather remote substituent R is questionable". They regrettably failed to notice in Table 1 and Fig. 1 of our manuscript that the C-Mg-C bond angle goes from 107.8°to almost collinear (179°) by the steric effect of the remote substituent.
Finally, we must also refer to Pan and Frenking most unjustied criticism: "[.] the reported bond dissociation energies (BDEs) given in Tables S5-S7 of the ESI of their work show very similar BDE values at different levels of theory, which question the value of the R (%) data". In fact, the R (%) values are obtained from the EOS analysis of the multireference CASSCF wavefunctions. There is no relationship between the BDEs and the R (%).
Pan and Frenking also confound some statements and data we provide in the introductory part of our original manuscript, which merely served as hints to illustrate that the original Be(0) assignment might not be appropriate. This background information was not implying that Be(0) assignment must be incorrect, but was rather highlighting the uncertainty in previous assessment that prompted the more comprehensive analysis we performed in our work.
From their comments it appears that we actually share the view on most of the fundamental aspects of the oxidation state concept, but again, they distorted our statements. First of all, we agree with Pan and Frenking that merely considering the atomic electronegativities of the contact atoms to decide the fate of the electrons of a bond is troublesome. Some of us have already discussed this issue in detail, 10 and that is precisely why we advocate for alternative approaches to simplistic IUPAC algorithms, such as EOS. Now, we did not question the existence of dative bonds between the ligands and Be/Mg on the basis of their respective electronegativities as Pan and Frenking insinuated. This is clearly not the point. The point is, in the context of IUPAC's ionic approximation (IA), any bond (dative or covalent) will have the corresponding electron pair assigned to either Be/Mg or C, and in the most straightforward application of IA, it should be assigned to the most electronegative atom, 20 which is the carbon atom. This is what stems from EOS analysis in the cases in question such as Be(cAAC Dip ) 2 , as illustrated in Fig. S37 of our original manuscript. The occupation of the s-type EFO on the ligands (0.85) is much higher than that of the s-type EFO on Be (0.18), which results in the electron pair of the s M-L bond unambiguously assigned to the ligands.
What is really relevant from the OS point of view is the assignment of the p electrons of the system. In the Be(0) picture, the electrons of the 3c-2e p-bond are assigned to the less electronegative atom, at odds with IUPAC's simple rules. 20 Since the latter are not foolproof we conducted EOS analysis, which in this case is crystal clear: as shown in Fig. S37 of ref. 2, the occupation of the p-type EFO on Be (0.12) is much smaller that of the p-type EFO of the ligands (0.44), which keep the electrons, resulting in a Be(+2) picture.
In our original manuscript, we also explicitly stated that Ponec and Cooper had also questioned the fate of the p-electrons using domain-averaged Fermi-hole (DAFH) analysis for a Be(cAAC H ) 2 model system. 21 According to Ponec and Cooper, the electron pair formally originally occupying the Be(2p p ) orbital contributes symmetrically to such an extent into the originally empty p orbitals on the carbene-like moieties that there is only a very small residual population on Be. In fact, the DAFH analysis reveals a contribution to the p bonding of 0.14 from Be and 0.95 from each of the ligands (see Fig. 1 of ref. 21 ). This is in full agreement with our results and with a Be(+2) formal valence state.
Pan and Frenking also appear to estimate the contribution of an AO to a given MO just by looking at the orbital plot. They focus on a p-type natural orbital with a natural occupation of 1.80. They claim without any proof that the largest coefficient of that orbital is at the Be atom. Then, according to IUPAC's IA, the electron pair should be assigned to Be. They quote a statement from ref. 20 that in our view represents a naïve and ambiguous point on the relationship between AO coefficients and OS. This issue has already been discussed elsewhere by some of us, 22 so we will only provide a quick answer to Pan and Frenking. First of all, the orbital they refer is a natural orbital corresponding to our CASSCF description of Be(NHC Dip ) 2 , not Be(cAAC Dip ) 2 . For the latter, the relevant orbital is given in Fig. S15 of our original manuscript, 2 with a natural occupation of 1.64. It is not so obvious which is the contribution of a given fragment holding a set of non-orthogonal AOs to a molecular or natural orbital. One may consider an orthogonal basis instead, such as a Löwdin's variant implemented in pySCF. 23 The orbital coefficients of the 2p z and 3p z AOs of Be are 0.399 and −0.096 for Be(NHC Dip ) 2 , and 0.291 and −0.095 for Be(cAAC Dip ) 2 . For comparison, just the coefficient of the 2p z AOs of each of the two contact carbon atoms of the ligands is 2 × 0.495 and 0.541 and 0.512 (nonsymmetric orbital), respectively. The population of these natural orbitals on Be (i.e. sum of the square of the coefficients on an orthogonal basis) is merely 0.18 and 0.11, respectively, not too different from the DAFH and EFO occupations discussed above. Looks can be deceiving.
We do concede to Pan and Frenking that the excitation energies of Mg that we discussed in the Introduction of ref. 2 are plainly wrong. That was a regrettable mistake and we are grateful to them for pointing it out, despite being irrelevant to the discussion. Pan and Frenking used the ionization potential and electron affinity of the fragments. They took the electron affinity of cAAC Dip from a PBE0/def2-TZVP calculation without Zero Point Vibration Energy (ZPVE) correction, 24 which does not meet the denition of molecular electron affinity 25 (it is not 18.4 kcal mol −1 but 14.5 kcal mol −1 ). Surely, the ionization energies required to achieve a real (not formal!) Be 2+ species are much higher than those needed to reorganize the valence electron in a Be(0) picture, but this does not preclude the possibility of a formal Be(+2) OS. Otherwise, how to justify the O(−2) formal OS in H 2 O, where the formation of O 2− (which is unstable, by the way) from O(0) is endothermic by +167.8 kcal mol −1 (ref. 26 and 27) and the ionization potential of H is +313.3 kcal mol −1 ? 28 Surely, Pan and Frenking are not challenging the formal OS assignments in the water molecule.
Another set of criticism from Pan and Frenking points towards our illustrative EDA calculations and EDA-NOCV analysis of these Be and Mg complexes, reported in the ESI of our original work. 2 Pan and Frenking stated: "we repeated the calculations at the same level of theory as the authors but with inclusion of dispersion interactions, which had not been considered. We found that the DE orb values using singly charged and neutral fragments become nearly equal, which means that Be(cAAC Dip ) may be described with the oxidation states Be(0) or Be(+1) but clearly not with Be(+2)." Unfortunately, these calculations were neither reported nor used in their comment, and instead they reused our original Table S8 as Table 1 in their comment, which in fact always included dispersion corrections. First of all, assuming they did not reoptimize the geometry, the inclusion of dispersion will only provide a new stabilizing term DE disp , with exactly the same DE orb values. Secondly, they complain about our notation, where the total orbital interaction term DE orb-corr is presented as a sum of DE orb and DE orb-HF contributions, and they erroneously assume that we used them for our interpretation. Taking into account that the analysis is performed at the B3LYP-D3(BJ)/ TZ2P level of theory, the latter corresponds to the contribution from 20% of exact-exchange of the density functional approximation. The reason the two terms are given explicitly is connected to the fact that the current implementation of EDA-NOCV does not decompose the non-local HF-like exchange contribution to DE orb . 29 Such information has been described in dedicated technical literature and also could be easily checked in the outputs of their own calculations, or even in former tables reported by Pan and Frenking. 30,31 Finally, and more importantly, Pan and Frenking claim that our analysis of the electron ow of the NOCVs is wrong just by how the deformation density looks like in our Fig. S41. They ignore the actual contributions of the symmetry-adapted fragment orbital (SFO) basis to the deformation density, which are also provided in our Fig. 1.
Let us repeat here our point concerning the p channel from the NOCV analysis. For the sake of simplicity, we will consider the total a and b contributions together. For each channel k, the deformation density is written as where n k is the eigenvalue (in absolute value) and the pair of NOCVs are expressed in the (orthogonal) SFO basis {c} by the coefficients c k . Integrating the deformation density leads to which readily affords a decomposition in terms of the SFO basis. The most signicant contributions of the SFOs can be found in the output of a typical NOCV run. When using a Be 0 reference, the corresponding eigenvalue is ca. j1.5j, meaning that an electron ow of 1.5 e is associated with this particular channel of the deformation density. In this case, it can be seen that the 2p z -type SFO on Be exhibits a very large and negative contribution of −1.49 e, meaning that electrons move from this orbital to those with positive contribution (ligands, in this case). This SFO contribution essentially fully accounts for the whole charge transfer of this channel. When using a Be 2+ reference, the eigenvalue of the deformation density is j1.8j, but the contribution of the 2p z -type SFO on Be is merely +0.42 e. The positive contributions from other SFOs of the Be moiety are negligible, so ca. necessarily +1.4 e move into the ligand (and from the ligands, as indicated by the ligand-centered larger negative SFO contributions). Regrettably, Pan and Frenking strongly relied on our selected isocontour value (our Fig. S41 of the original manuscript and Fig. 2 in their comment) although the orbital interaction shows a Be contribution only at low isocontour values, while for bigger isocontour values it shows only the participation of the carbene carbon of the cAAC ligands (Fig. 2). This intrafragment charge-transfer process can be even more clearly seen with the NOCV orbitals, where the major contribution comes from the p-orbital of the cAAC ligands (Fig.  1). A similar analysis can now be carried out not in the orbital space but in the real space, using a Hirshfeld-type partitioning for the entire fragments (A), 32 which affords the same picture albeit with a much smaller charge transfer (−0.78 e on Be for the Be 0 reference and +0.04 e on Be for the Be 2+ reference) (Fig. 1). Hence, the actual numbers, and not a mere picture, evidence that there is much less electron ow between fragments using the Be 2+ reference, thus invalidating the Pan and Frenking claim. For completeness, we have performed the EDA-NOCV analysis with the correct wavefunction also for the rather odd Be 1+ fragmentation (see Table 1). The NOCVs of the p-bonding channel using a Be 1+ fragmentation indicate an intermediate situation between using Be 2+ and Be 0 , as might be anticipated. The eigenvalue of the deformation density is j1.54j, while the contribution of the 2p z -type SFO on Be is −0.51 e, slightly larger than the value using Be 2+ discussed above (+0.42 e). In the Be 1+ case, however, it originates from a partial cancelation of a (−0.75 e) and b (+0.24 e) contributions, thus indicating a much signicant charge-reorganization than when using Be 2+ as reference. Yet, from an energy perspective, Be 0 and Be 1+ show lower DE orb values than the Be 2+ fragmentation scheme. Does this necessarily mean that the oxidation state of Be cannot be Be(+2), while all other evidence clearly points toward this situation? Has the criterion of the lowest orbital interaction term been carefully tested in the past? In this sense, some of us considered the simplest case of H 2 O when evaluating an OS assignment method based on the position of the centroids of localized orbitals with respect to the atoms. 33 We found that the method failed to describe H 2 O as formally H(+1) and O(−2). As will become clear later on, this kind of proof of concept calculations were regrettably skipped when putting forward EDA's orbital interaction energy as a general criterion to establish appropriate valence states and OS (vide infra).
Aer having addressed Pan and Frenking's inaccurate statements, let us go back to square one. As mentioned above, the closed-shell KS-DFT wavefunction of many of these systems are either not stable or a lower in energy broken-symmetry solution is found. In the case of Be(cAAC Dip ) 2 , the corresponding hS 2 i value is 0.57, which cannot be ignored at all. Indeed, this is supported by the computed D1 diagnostic, 34 which exceeds in all the cases the reference value. Consistently, the weight of the closed-shell conguration is for all the complexes below 0.85. Local spin analysis 35 thus reveals a signicant diradicaloid character, with local spins centered at both cAAC Dip ligands. That is why we claim that these systems have diradicaloid character, not because they exhibit a small singlet-triplet gap, as Pan and Frenking stated. The CASSCF description also conrms this picture, both from the signicant occupation of the antibonding natural orbital of the Be(cAAC Dip ) 2 p bond (0.36) and from the local spin analysis values reported in Table 1 of our original paper. 2 This situation is quite similar to that of the largely discussed bonding in the [NaBH 3 ] − cluster. 36 In our earlier work on that subject, 36 we already showed that "EDA cannot distinguish an electron-sharing interaction from a spin-polarized one (diradicaloid)", at least by only looking at the energy components. Pan and Frenking do not confront the statement and clearly underestimate the implications of this fact, by stating that these are subtle differences and fall outside the scope of EDA. We could not agree more with their own statement: "the EDA is as good as the quantum theoretical method on which it is based". 1 We in fact extend it to any wavefunction analysis method. That is precisely why one should make sure that an appropriate wavefunction is considered. At least one can check whether it has internal instabilities or not. But how can this not be relevant for OS assignment from wavefunction analysis? If a given bond between A and B is "electron-sharing", the electron pair will be allocated to either A or B depending on the direction of the ionic approximation. If there is enough spin polarization of the bond, a more reasonable result might be the homolytic distribution of electrons, 37 since in the extreme case (a pure diradical) there is no "sharing" of electrons. The structure of the singledeterminant closed-shell wavefunction and that of a brokensymmetry one can be very different, and this can even be taken into consideration for the EDA analysis (Fig. 3). In contrast, they justify their systematic misuse of EDA with "[.] it is odd to criticize the method for not providing information about a particular property that is not the target of EDA." Should we stop caring about using the correct wavefunction?
Notably, Pan and Frenking brought the case of the Ca-CO 2 complex in their comment, 39 as an example where EDA-NOCV works for a system with open-shell singlet ground state character. Indeed, this is not their best example since one can clearly notice in their original paper that they have in fact used once more a closed-shell unstable wavefunction in their EDA-NOCV analysis, ignoring their own described diradical character. 39 We have carried out the correct calculations, and we also present herein the obtained results. Table 2 gathers the outcome of the EDA for the closed-shell (unstable) and openshell wavefunctions for the situation where the fragments are a Ca + (4s 1 ) cation in the doublet ground state and a doublet The value in parentheses gives the percentage contribution to the total attractive interactions DE elstat + DE orb-corr + DE disp . b The values in parentheses give the percentage contribution to the total orbital interaction DE orb-corr .    with R (%) = 93.6, which is the assignation Frenking reported by analyzing the wrong wavefunction. It is evident that Pan and Frenking are only stacked in the criterion of the lowest orbital interaction term, which looks more like a "credo" than a carefully tested indication of the best electronic reference situation. It is a good exercise to track back where their recited creed "those fragments which give the lowest energy change for the orbital interactions are the most suitable species to describe the bonding interaction" comes from. According to them, this is a proven fact; however, their citations only go to articles with the same phrase and no proof. To the best of our knowledge, the rst mention of this was given in 2007, 40,41 but there is no sign of a thorough systematic study of this matter, only a belief. [42][43][44][45][46] One should not forget the memorable Feynman's quote: one must not verify an idea using the same data that suggested the idea in the rst place. It is interesting to point out that some authors have already questioned its validity. 6,47 Also, as repeatedly stated by the original developers of the EDA, the orbital interaction term includes both charge-transfer and polarization contributions, also within the fragments as shown here, which can certainly blur the analysis. 48-51 Some of us have been diagnosing the problems associated with such an energetic criterion when it comes to OS assignment in general (notice that it is never considered in any of the IUPAC technical reports), and found that EDA-NOCV can still be a valuable tool if the interfragment electron ow is considered instead. It is out of the scope of the present response to discuss the details, which will be given elsewhere soon, 52 but we can already disclose a quite surprising result. There is no need to use the exotic Be(cAAC Dip ) 2 system to observe that the orbital interaction is not a good indicator of the best electronic reference. As shown on Table 3, BeH 2 would be in fact best described as Be(0) and H(0) by the lowest orbital interaction energy criterion. Are Pan and Frenking also challenging the oxidation state of BeH 2 ?

Conclusions
Our assignment of Be(+2) and Mg(+2) oxidation states in Be(cAAC Dip ) 2 and Mg(cAAC Dip ) 2 is fully consistent with our data. Although Pan and Frenking have claimed a critical inspection of our data, their arguments are based on supercial visual inspection of our gures rather than a thorough numerical analysis. Additionally, we showcase their recent awed application of EDA to diradicaloid species Ca-CO 2 . We provide evidence that the criterion of the lowest orbital interaction energy using EDA to determine the valence state of fragments is intrinsically unreliable, since it already fails for a simple system such as BeH 2 . As the oxidation state is a density based property, NOCV analysis can be a more reliable indicator of the best reference state. EDA-NOCV is without doubt a very useful method for bonding analysis, yet it should be used with appropriate care.