Investigating István Mayer's “ improved ” definitions of bond orders and free valence for correlated singlet-state wave functions

Chosen from István Mayer's very impressive canon of important work on bond orders and related quantities, we explore a paper from 2012 which introduced a “ pseudo spin density matrix ” for correlated singlet-state wave functions, leading to “ improved ” definitions of bond orders and free valence. Examining such “ improved ” bond orders for the singlet ground states of H 2 , N 2 and CH + we find that all of them exhibit sensible geometry dependences. Mayer's free valence index works well for H 2 and N 2 ; the asymptotic behavior for CH + turns out to be slightly more compli-cated but can easily be explained. Using B 2 H 6 to examine three-center bonding, a simple generalization of Mayer's approach produces numerical results close to those based on the “ pseudo spin density matrix. ” As expected, the various “ correction ” terms remain small, albeit some of them are larger for the multicenter indices of ground and excited singlet states of benzene, S 2 N 2 and square (D 4 h ) cyclobutadiene.


| INTRODUCTION
As is well known, a central theme of István Mayer's distinguished scientific career was the extraction of useful chemical information from molecular wave functions. Among his principal requirements for such methodologies were that they should be well-defined and produce results that make proper chemical sense. He was of course particularly active in the development of successful measures of bond orders and related quantities. He published an interesting personal account of his work up to 2006 on the development of bond orders [1], with various subsequent studies, as well as work on energy components, being included in his 2017 book [2]. He identified, in the intervening years, a specific issue with his usual definition of bond orders for correlated singlet-state wave functions, publishing in 2012 his "improved" definitions of bond orders and free valence for correlated singlet-state wave functions [3].
The well-chosen specific example in that 2012 paper [3] and subsequently in his 2017 book (on pages 99-100) [2] is the dissociation of the singlet ground state of ethene into two triplet methylene fragments. Whereas the calculation of bond orders using the two open-shell CH 2 subsystems does of course involve the spin density associated with the unpaired electrons, the overall singlet state of H 2 C CH 2 has zero spin density and so there can be no such net contribution to the bond orders. The effect on the calculated values turned out to be marginal in this case but, characteristically, Mayer did of course note the conceptual inadequacy of being able to obtain different results for two triplets as opposed to one overall singlet [3]. Such issues do of course potentially arise whenever an overall singlet state dissociates into fragments with non-zero spin.
Mayer's solution to this general problem required the introduction for correlated singlet-state wave functions of what he termed a "pseudo spin density matrix" R, leading to so-called "improved" definitions of bond orders [3].
A key purpose of the present work is to explore further István Mayer's "improved" definitions of bond orders [3] as well as his related Rbased definition of free valence [2,3]. As diatomic examples we consider first the breaking of the formal single and triple bonds in H 2 and N 2 before examining the corresponding process for a heteronuclear example, namely the CH + ion. We then move on to three-center bonding using as an example the familiar case of B 2 H 6 , for which we compare results obtained with the matrix R with those from an alternative scheme. Along the way, we also examine the characteristics for these various molecular systems of Mayer's proposed R-based definition of free valence; this leads to a mix of predictable and unexpected outcomes. Finally, we employ ground and excited singlet states of benzene, S 2 N 2 and square (D 4h ) cyclobutadiene to investigate analogous "improved" definitions of multicenter indices for such aromatic and antiaromatic systems.

| THEORETICAL BACKGROUND
For certain special cases, such as a single determinant closed-shell restricted Hartree-Fock (RHF) N-electron wave function, the product DS of the usual charge density matrix D and the corresponding overlap matrix S is duodempotent, that is, (DS) 2 =2DS. It thus follows for k =1,2,3,… that: Focusing first of all on Wiberg-Mayer bond orders, we may write: in which P σ is the charge density matrix for σ spin, with D = P α + P β , and the notation μA signifies that the relevant summation is restricted to the basis functions associated with atom A.
Making use of the definition of the spin density matrix, P S = P α − P β , we may trivially rewrite Equation (2) in the following equivalent form: In practice, of course, DS is rarely duodempotent for more general wavefunctions, but we may still use essentially the same partitioning of tr[(DS) k ] to define Wiberg-Mayer and higher-order bond indices.
A fundamental problem with using Equation (3) for the dissociation of ethene (as mentioned in the Introduction) arises because the spin density matrix P S is necessarily null for correlated singlet-state wave functions. István Mayer's solution to this problem required the introduction of a "pseudo spin density matrix" R, to be used instead of P S [2,3]. Mayer's definition of R can be traced back to the introduction by Takatsuka et al of a quantity which may be termed the distribution of effectively unpaired electrons [4,5]. Generalizing such definitions, we may define a hierarchy of matrices Q (k) according to the value of k =1,2,3,… such that: In practice, Q (k) can conveniently be generated as follows: Not only does the resulting matrix Q (k) for k = 2 coincide with Mayer's matrix R, but our scheme for generating it also matches the procedure that he presented [3]. Mayer's "improved" bond order, which is obtained by substituting R for P S in Equation (3), then takes the form [2,3]: Not only does this new definition solve the problem that Mayer had identified for the dissociation of ethene [2,3] but it has the added advantage that the sum of all possibleW R AB values for a system with N electrons is now necessarily equal to 2N, as in the case of a single determinant closedshell RHF wave function.
The same matrix R also performs a role analogous to the spin density matrix in Mayer's related definition of a free valence index for correlated singlet-state wave functions. This quantity,F A , which takes the form [2,3]: is designed to be a measure of the extent to which the electronic structure around a given atom deviates in the correlated case from a situation with only doubly occupied natural orbitals. (Part of István Mayer's justification [2,3] for Equation (7) involves an expression [6,7] for the total valence, V A ; we note that this definition of V A appears in his 2017 book [2] as equation (6.57) but that it was unfortunately misprinted in the 2012 paper [3], with some symbols missing from his equation (14).) The Mulliken-like (Hilbert space) partitioning that is implicit in many of the above equations does of course restrict the types of basis set that can reasonably be used. In order to produce results which make proper chemical sense, it is typically important to use basis sets with a high degree of atomic character, as Mayer explains in his 2017 book (on page 14) [2]. Among others he identified the 6-31G** basis set as being suitable. On the other hand, contemporary large basis sets, typically including diffuse functions, are less likely to be well-suited and this does of course also imply that there is no obvious basis set limit for quantities obtained with this type of approach.
As is well described in chapter 9 of István Mayer's 2017 book [2], an appealing alternative which circumvents such basis set issues is to perform the analysis in three-dimensional space with suitably defined atomic domains, whether overlapping ("fuzzy") or disjunct. An obvious example of the latter is the widespread use of quantum theory of atoms in molecules [8] (QTAIM) domains, Ω A . This is the approach that we have adopted for almost all the results described here.
When working in the basis of (real) orthonormal natural orbitals ϕ K with occupation numbers n K , the matrix D is of course diagonal, with elements n K , and S is a unit matrix. The matrix Q (k) also turns out to be diagonal, with elements: and we note that such values are clearly zero when n J is 0 or 2.
The QTAIM-generalized [9] "improved" Wiberg-Mayer index then takes the following relatively simple form [10,11]: in which KL ji ΩA D denotes a domain-condensed overlap integral: The free valence index,F A , as defined in Equation (7), also adopts a simple QTAIM-generalized form, namely: Given that R JJ =[n J (2 − n J )] ½ , it is of course now rather easy to see from Equation (11) that this index does indeed quantify deviations from double occupancy.
We turn now to the definitions of the three-center indices which can be obtained from the partitioning of 1 /4 tr((DS) 3 ). Starting with the relevant QTAIM-generalized expressions for the case of non-zero spin and then recognizing again that the spin density identically vanishes everywhere for a correlated singlet-state wave function, it does of course make sense for such singlet states to replace all occurrences of P S with the matrix R. As has been shown before [11,12], this strategy leads to the following expression for an "improved" three-center quantity which we denote X R ABC : However, we notice in this case that [11]: so that there is no direct relationship of the sum of all possible X R ABC values to the total number of electrons, N. This potential drawback has prompted us to propose the following admittedly ad hoc generalization of Equation (9) to the three-center case [10,11]: It is straightforward to show that the sum of all possible X Q ABC values is equal to N, essentially by construction. Furthermore, this very simple expression for the three-center quantities X Q ABC can of course easily be generalized to higher numbers of centers. One could even trivially generalize the definition of F A in Equation (11) to k > 2, again using the appropriate Q (k) matrix, but it is not clear whether the resulting quantities would be of any practical use.

| COMPUTATIONAL DETAILS
We employed full configuration interaction wave functions for H 2 and full-valence complete active space self-consistent field (CASSCF) calculations for both of N 2 and CH + to study the variation with nuclear separation, R, of the Wiberg-Mayer bond order and the free valence index F A for each system. All the required wave functions were calculated in D 2h or C 2v symmetry, as appropriate, using standard cc-pVTZ basis sets. In practice, the wave functions for H 2 were generated by means of "2 electrons in 28 orbitals" CASSCF calculations, those for N 2 involve "10 electrons in 8 orbitals" [11,13] and those for CH + (X 1 Σ + ) involve a "5 electrons in 4 orbitals" construction. For each of these three molecular species, CASSCF calculations were carried out for a range of values of R using the MOLPRO [14][15][16] quantum chemistry package. We do not expect expansions of the CASSCF active spaces or extensions of the basis sets to lead to significant differences in any of the R dependences of the quantities we examine here.
For investigating three-center indices we used the same wave function for B 2 H 6 as in previous work [11,13], namely full-valence CASSCF ("12 electrons in 14 orbitals") with a standard cc-pVTZ basis set. An obvious advantage of using the same natural orbital expansion as in a previous study is that this enables direct comparisons of results. For much the same reason, when investigating multicenter indices (k > 3), the wave functions for ground and excited singlet states of benzene, S 2 N 2 and square (D 4h ) cyclobutadiene were taken from a recent study of these various aromatic and antiaromatic ring systems [17]. This means in practice that the wave functions were calculated at the CASSCF(6,6) and CASSCF(4,4) levels for C 6 H 6 and square (D 4h )C 4 H 4 , respectively, using the 6-311++G(2d,2p) basis set, and at the CASSCF(22,16) level for S 2 N 2 , using instead the cc-pVTZ basis set. A few wave functions were also generated for these ring systems with Gaussian 09 [18] using the 6-31G** basis set so that we could use Mulliken-like (i.e., Hilbert space) partitioning for comparison.
The various bond indices and free valences reported in the present study were calculated using our own programs. Domain-condensed overlaps required for the QTAIM-generalized approach were generated using the AIMAll program [19].

| RESULTS AND DISCUSSION
We compare in Figure 1A the dependence on R for H 2 of the QTAIM-generalized "original" and "improved" definitions of the Wiberg-Mayer bond order. For R = 0.75 Å, near R e , the "improved" value is W R HH 0 =0:961, with the R-based correction term in Equation (9) accounting for only 0.012. Indeed, the curves obtained with the two definitions remain very close to one another over the full range of nuclear separations that we considered (0.4-20 Å). Both values, and thus also the R-based correction term, do of course go to zero at large R. The monotonic decrease of the Wiberg-Mayer bond order with increasing R can be seen ( Figure 1A) to feature a point of inflection near 1.84 Å. Looking instead at the free valence index, the values increase monotonically from F H = 0.059 at 0.75 Å toward unity at large R (see Figure 1A).
The general form of the observed R dependence of W R HH 0 and of F H does of course meet our chemical expectations, with the two curves being almost mirror images of one another, so that the curves for W R HH 0 and 1 − F H run essentially parallel to one another. This does though have the obvious consequence in the present case that an examination of the R dependence of F H does not provide new information.
The corresponding results for the dissociation of N 2 into two N( 4 S) atoms, as depicted in Figure 1B, tell a somewhat similar story. For R = 1.1 Å, near R e , the "improved" Wiberg-Mayer bond order is W R NN 0 =2:819, with the R-based correction term in Equation (9) accounting for only 0.027. As was the case for H 2 , the values decrease monotonically to zero, with the two curves remaining close to one another. There is again a single point of inflection, this time near 1.81 Å, whereas we might have hoped to see a first one for the breaking of the double π bond and then a second one at larger R for breaking the σ bond.
We find that the corresponding value of F N increases monotonically from 0.265 at 1.1 Å toward the expected value of 3 at large R, with the curve for F N being almost the mirror image of the one for W R NN 0 ( Figure 1B). Indeed, the curves for W R NN 0 and 3 − F N run almost parallel to one another so that, just as was the case for H 2 , the R dependence of F N does not provide new information.
As can be seen from Figure 2A, the R dependence of the Wiberg-Mayer bond order for the formal single bond in the 1 Σ + ground state of CH + follows the same basic pattern as for H 2 and N 2 . The "improved" Wiberg-Mayer bond order for R = 1.1 Å, near R e ,isW R CH =0:914 with the R-based correction term accounting for only 0.009. The correction again turns out to be small over the full range of R that we have considered. It does though turn out that the corresponding values of F A for the carbon and hydrogen centers do not follow so closely the behavior that was observed for H 2 and N 2 .
Given that the first ionization energy of the C atom is somewhat lower than that of H, the dissociation limit in cylindrical symmetry is C + (2s 2 2p z ) + H(1s). Such a situation can be realized by natural orbital occupation numbers of 2, for C + (2s 2 ), and then two values of 1, for C + (2p z ) and H(1s), and zero occupancies for both of C + (2p x ) and C + (2p y ). As such we could reasonably expect the limiting values of F A for either of the atomic domains to be unity. As is shown in Figure 2A, this is indeed what is observed for the hydrogen center, for which the value of F H increases monotonically from 0.030 at 1.1 Å toward unity at large R. The corresponding value for the carbon center at R = 1.1 Å is 0.537 and we again observe a monotonic increase with increasing R, but clearly toward values which significantly exceed unity (see Figure 2A).  We constructed a simple model for this full-valence CASSCF description (see Figure 2B) by recognizing that the value of F C is given by the expression 2x(2 − 2x)+1+2x(2 − x), so that F C exceeds unity by 8x − 6x 2 . Clearly there will be a maximum value of F C = 11 3 11 / 3 when x = 2 /3, but this would of course be a somewhat extreme situation. Even so, the function 8x − 6x 2 does vary fairly rapidly even for small values of x.
Given that the actual values of x for CH + from the full-valence CASSCF calculations at large R is on the order of 0.05, our simple model predicts F C 1.385. Such values, significantly exceeding unity, are indeed what we have found. Far more straightforward behavior would undoubtedly be observed when using instead local spin analysis [20], for which the asymptotic values of hŜ 2 i A for the C and H centers would both be 3 /4. Even though it turns out to be entirely straightforward to understand how the observed F C values on the order of 1.4 occur for formally C + (2s 2 2p z ), they certainly diminish somewhat our confidence in using such quantities in comparisons of the bonding situations in different molecules. Now, using B 2 H 6 as our example, we consider also the three-center analogs of Wiberg-Mayer bond orders which correspond instead to partitioning of 1 /4 tr[(DS) 3 ]. Before doing so, we look first at the hydrogen atoms in this molecule, for which we find F A values for the bridging (H b ) and terminal (H t ) domains of 0.061 and 0.044, respectively. We notice that these small values are comparable to the value of 0.059 that we found for H 2 near its equilibrium geometry.
Instead of examining in turn the Wiberg-Mayer bond orders for each different pair of individual atoms as well as the corresponding threecenter indices for all possible triads, it proves more convenient to assemble the various B and H t QTAIM domains into two BH 2 domains. Accordingly, we report in Table 1 our free valence indices, Wiberg-Mayer bond orders and three-center indices for the two BH 2 and two H b domains.
Given the invariance of any of the various X ABC quantities to permutations of the indices it is convenient, when all of three indices are different, to quote in each case a (single) value of X(A,B,C)=3!X ABC , as we have done in Table 1. In keeping with our traditional notions of three-center two-electron bonding across the bridging hydrogen atoms, we observe from Table 1 that the largest three-center index occurs for the BH 2 ,H b , B 0 H 2 triad. There is though also a nontrivial value for the BH 2 ,H b ,H b 0 triad. A further observation from these various values is that the correction terms involving the matrices R and Q (3) are small, with it making hardly any difference whether we use values of X R (A,B,C)orX Q (A,B,C).
The largest of the "improved" Wiberg-Mayer bond orders, W R AB =0:502 (see Table 1), occurs for the BH 2 ,H b pair, as might have been expected, with the R-based correction term accounting for just 0.013. We do though also observe smaller Wiberg-Mayer bond orders, on the order of 0.2, for the BH 2 ,B 0 H 2 and H b ,H b 0 pairs. It is important not to suppose that these values necessarily indicate a degree of direct chemical bonding. As was noted, for example, by Ponec and Uhlik [21], the existence of genuine three-center indices for the ABC triad demands nontrivial indices not only for the AB and BC pairs but also for the AC pair. In his 2017 book [2], István Mayer described (on page 78) such an AC bond order in the absence of direct chemical bonding as resulting from a sort of interference, potentially enhancing molecular stability via an attractive exchange contribution.
As our final group of examples, we now use correlated wave functions for ground and excited singlet states of benzene, S 2 N 2 and square (D 4h ) cyclobutadiene to investigate "improved" definitions of multicenter indices (k > 3) for such aromatic and antiaromatic systems. Instead of considering the C 4 H 4 and C 6 H 6 rings in terms of separate C and H atoms, as was done in our recent study [17], we have chosen here to perform the analysis in terms of merged CH domains. We find that the values of the free valence index F CH for these hydrocarbon rings all lie in the range 2.81-2.98, with no sign of any link to whether a given state is aromatic or antiaromatic.
Our various four-and six-center indices, with and without Q-based corrections, are reported in Table 2. Note that these values include the terms for all permutations of the indices and, as in our previous study [17], they have been calculated using spin-orbitals instead of orbitals, enabling direct comparisons with the previous work. The resulting QTAIM-generalized quantities for C 6 H 6 turn out to be slightly larger for six CH units than for just the six C atoms [17], with the opposite being true for the corresponding states of square (D 4h )C 4 H 4 . These small differences do not modify the patterns of values for the various singlet states of either of these molecules.
For the sake of comparison, we also calculated multicenter indices and Q-based correction terms for some of these states using Mulliken-like (i.e., Hilbert space) partitioning with the 6-31G** basis set. In particular, we also report in Table 2 our Q-corrected values,X Q 6 , for the six C atoms in the S 0 ,S 1 and S 2 states of benzene.
We observe from Table 2, especially for C 6 H 6 and S 2 N 2 ,agreaterimpactoftheQ-based correction terms on the QTAIM-generalized indices than was the case for the diatomic molecules and for B 2 H 6 . Even so, the basic patterns of values for the various singlet states of each molecule do still turn out to be much the same as we have reported previously [17] so that the inability of these multicenter indices to distinguish properly between the aromatic and antiaromatic singlet states of a given molecule has not been remedied by the inclusion of the correction terms. Instead, we still observe the lack of any clear distinction between, for example, the six-center indices for the antiaromatic S 1 and aromatic S 2 states of benzene [17].
The Q-based correction terms for the six C atoms of benzene from Mulliken-like (Hilbert space) partitioning turn out to be somewhat smaller than are those from the QTAIM-generalized approach. Indeed, even when expressed to four decimal places, it is only for the S 2 state that they make any difference, increasing the overall value from 0.0037 to 0.0039. Clearly the values ofX Q 6 ( Table 2) again do not properly distinguish the S 1 and S 2 states of this molecule, so that our key conclusion is the same as before.
The underlying problem remains after the inclusion of these correction terms, namely the reliance on one-electron density matrices (and thus occupation numbers). Such density matrices do of course reflect the multi-reference character but they can, as has been shown previously [17], nonetheless obscure important details that are apparent in the CASSCF wave functions, such as the degree of partial diradical character.

| CONCLUSIONS
We have explored in the present work a paper from 2012 [3], chosen from István Mayer's very impressive canon of important work on bond orders and related quantities. It is clear that Mayer was not at all content that the numerical consequences of defects in some of his original definitions were marginal for the dissociation of the singlet ground state of ethene into two triplet methylene fragments [2,3]: there was nonetheless a conceptual inadequacy to be resolved. His solution required the introduction for correlated singlet-state wave functions of a "pseudo spin density matrix" R, which led to "improved" definitions of bond orders and free valence [2,3].
To explore these new definitions, we examined the breaking of the bonds in the singlet ground states of H 2 ,N 2 and CH + .The"improved" bond orders do show sensible geometry dependence even if, as is to be expected, the R-based correction terms are always small. Mayer's free valence index works well for H 2 and N 2 ; the asymptotic behavior in the case of CH + turns out to be slightly more complicated, but it can nonetheless easily be explained.
Treating B 2 H 6 as two BH 2 units linked by two bridging H atoms, we explored the behavior of three-center indices in this familiar example which transcends the usual Lewis model through the incorporation of three-center two-electron bonding units. A simple generalization of Mayer's approach, using a matrix which we have denoted Q (3) , was found to produce numerical results close to those from the R-based scheme. As is to be expected, all of these various correction terms remained small.
Moving on to four-and six-center indices, we examined values for ground and excited singlet states of C 6 H 6 ,S 2 N 2 and square (D 4h )C 4 H 4 .
The Q-based correction terms turned out to be larger than for the diatomic molecules and for B 2 H 6 . Nevertheless, they did not impact the basic inability of these multicenter indices to distinguish properly between the aromatic and antiaromatic singlet states of a given molecule.
There are many who might ignore marginal effects on calculated values but, whether it was in the 2012 paper [3] or indeed in many of his others, science benefitted greatly from István Mayer's dedication to developing well-defined and conceptually valid approaches which lead to results that make proper chemical sense.