Electronic Structure of Cubane-Like Vanadium–Nitrogen Cationic Clusters [V₄N₄]⁺ and [V₆N₆]⁺

Abstract

Density Functional Theory and Complete Active Space Self-Consistent Field (CASSCF) methodologies are used to explore the electronic structure of the cationic V–N clusters, [V₄N₄]⁺ and [V₆N₆]⁺, that have been identified in recent mass spectrometric experiments. Our calculations indicate that both clusters are based on cubane-like fragments of the rock-salt lattice. In the smaller [V₄N₄]⁺ cluster, the V–V bonding is delocalized over the tetrahedron, with net bond orders of 1/3 per V–V bond. In [V₆N₆]⁺, in contrast, the V–V bonding is strongly localized in the central V₂N₂ unit, which has a short V=V double bond. CASSCF calculations reveal that both localized and delocalized V–V bonds are highly multi-configurational.

1. Introduction

Small clusters of transition metal nitrides, M_xN_y, represent the simplest models for the active sites in the bulk nitrides that have an important role in heterogeneous catalysis [1,2,3]. Amongst these, vanadium nitride is a particularly important example, and catalytic applications of both the bulk material and thin films have been reported [4,5,6,7]. Metal–metal bonding remains very strong in bulk VN, to the extent that it is an extremely good conductor [8], and the solid also features a high concentration of nitrogen vacancies, giving a range of stoichiometries VN_x, where x ≤ 1 [9,10]. Indeed, low-valent binary compounds of the early transition metal elements are, in general, stabilized by extensive metal–metal bonding; examples include TiN [11] and VO [12], as well as the low-temperature form of VO₂ [13]. The nature of the surface sites on non-stoichiometric VN remains a matter of debate, but it is clear that V–V bonding has the potential to drive extensive surface reconstruction that may impact the catalytic performance [1]. In light of this extreme non-stoichiometry in the solid state, small molecular clusters with V:N ratios close to 1 (i.e., V_xN_y, x ≈ y) provide an important opportunity to explore the interplay of metal–metal bonding, structure, magnetism and catalytic competence of the vanadium nitride phase. In recent work by Hirabayashi and Ichihashi, particular members of the [V_xN_y]⁺ family were shown to be highly active catalysts for ammonia decomposition [14]. Clusters with V:N ratios close to 1.0 were formed by sputtering metal targets with accelerated Xe⁺ beams in the presence of N₂. Prominent peaks in the mass spectrum were assigned to [V₃N₃]⁺, [V₄N₄]⁺, [V₅N₅]⁺ and [V₆N₆]⁺, as well as the vanadium-rich species [V₃N₂]⁺, [V₄N₃]⁺, [V₅N₄]⁺ and [V₆N₅]⁺. The subsequent reactivity of these ions with NH₃ depends critically on their composition, and the equiatomic clusters [V₄N₄]⁺ and [V₆N₆]⁺ appear to be particularly inert, binding one or two molecules of NH₃ without showing any N–H bond cleavage. All other ions, in contrast, lose an equivalent of H₂ from the bound NH₃ to form [V_xN_y(NH)]⁺. The relative inertness of [V₄N₄]⁺ and [V₆N₆]⁺ led Hirabayashi and Ichihashi to propose structures based on cubane-like architectures, which represent fragments of the rock-salt lattice of the parent nitride—VN. Cubane-like motifs are, in fact, relatively common in gas-phase clusters of binary compounds of the alkali and alkaline earth metals and can be viewed as intermediate stages in the crystal growth process [15,16,17]. (NaCl)_x clusters, for example, have been shown to adopt rock-salt motifs [18], as have a range of alkaline earth metal-oxide clusters, such as [M_xO_x]⁺, [M_x+1O_x]⁺ and [M_xO_x]²⁺ (M = Ca, Mg). Amongst the transition metal elements, the nitrides [Ti_xN_x]^0/+ (x = 4, 6 and 9) show similar structural features [11,19,20].

Our interest in the electronic structure of these cubane clusters lies in the extent to which they are stabilized by metal–metal bonding. Stable, coordinatively saturated cubane clusters are known across much of the transition series (see Figure 1 for selected examples), and many have important roles in biology. Fe₄S₄ clusters, for example, are ubiquitous in biological electron transfer [21] and were amongst the first to be studied using the broken-symmetry approach pioneered by Noodleman [22]. More recently, the remarkable structure of the Fe₇Mo double cubane in the FeMo cofactor in nitrogenase has been revealed [23,24], as has the importance of Mn₃Ca cubane motifs in water oxidation [25]. The presence of multiple bridging ligands between the metals in all these clusters means that whilst metal–metal bonding is possible, it is not absolutely essential for stability, and as a result, the bonding in cubanes spans the entire spectrum from strong covalent to weak exchange coupling. The cyclopentadienyl-capped clusters Cp₄M₄E₄ (M = Cr, Mo; E = O, S) [26] are a case in point wherein delocalized covalent bonding in the molybdenum systems gives way to weak exchange coupling in the chromium analogues.

In this paper, we explore the electronic structures of two of the clusters identified in the mass-spectrometric experiments of Hirabayashi and Ichihashi, [V₄N₄]⁺ and [V₆N₆]⁺, both of which have been proposed to adopt cubane-like architectures on the basis of their reactivity (or lack thereof) towards ammonia dehydrogenation. The coordinatively unsaturated nature of the metal centres leads to a relatively small splitting within the d-orbital manifold, and hence to a very rich electronic structure with many accessible electronic states. From a computational perspective, the description of metal–metal bonding therefore represents a very considerable challenge, and the extent to which single-determinant approaches, such as density functional theory, are appropriate remains to be established. Gorelsky [27] has explored this issue in the context of a vanadium dimer, V₂(C₅H₅)₂(pentalene), where both the geometry and the singlet–triplet splitting were found to be well described by the gradient-corrected (GGA) BP86 functional and also by the meta-GGA TPSS functional. The inclusion of exact Hartree–Fock exchange in the B3LYP functional, for example, causes a very substantial weakening of the V–V bond, which appears inconsistent with the available experimental data. In recent work, we identified a number of other cases where the inclusion of exact exchange in the functional provides a description of metal–metal bonding that is qualitatively different from the picture derived from GGA or meta-GGA functionals [28]. Here, we report results obtained using a range of different density functionals and compare these to the electronic structure picture that emerges from a parallel series of calculations performed with the CASSCF methodology.

2. Results

2.1. Electronic Structure of [V₄N₄]⁺

For [V₄N₄]⁺, we identified two closely-spaced states, both of which are based on a tetrahedral architecture (Figure 2). The first of these is a ⁴A₂ state with rigorous T_d symmetry: the six V–V bond lengths are identical at 2.65 Å (using the TPSS functional [29]). In the alternative sextet state, ⁶A₁, the tetrahedron is slightly elongated along one of the two-fold axes, giving four long (1.93 Å) and eight short (1.85 Å) V–N bonds and overall D_2d symmetry. The average V–V bond lengths in the ⁶A₁ state are somewhat smaller than those in ⁴A₂ (2.60 Å vs. 2.65 Å), while the average V–N bond lengths are marginally longer (1.88 Å vs. 1.86 Å). The origin of these structural patterns is apparent in the Kohn–Sham orbital diagram, also shown in Figure 2. In the perfect tetrahedron, the V $d_{z^2}$ orbitals form a basis for a₁ and t₂ representations (the z axis is defined locally here as the vector passing through each V atom and the centre of mass of the cluster). The a₁ combination is bonding with respect to the V–V interactions but also strongly V–N antibonding, and as a result, it lies above the weakly V–V bonding t₂ set. The $d_{xy}/d_{x^2-y^2}$ orbitals, in contrast, are aligned tangentially around the surface of the tetrahedron, and form a basis for e + t₁ + t₂ representations, the doubly degenerate e orbital being strongly V–V bonding. In the ⁴A₂ state, the two components of the 4e orbital are doubly occupied, while the three components of 12t₂ each contain a single electron, giving rise to a perfectly symmetric structure with six equivalent V–V bond lengths of 2.65 Å. Promotion of a single electron from the 4e orbital into 9a₁ results in an orbitally degenerate ⁶E state, and hence to a first-order Jahn–Teller instability that drives the distortion to D_2d symmetry (and thus to the ⁶A₁ state). The population of the 9a₁ orbital also strengthens the V–V bonds while at the same time weakening the V–N bonds, leading to a global contraction of the V₄ tetrahedron.

The structural trends identified in the previous paragraph are broadly independent of the choice of functional (

Table 1. Relative energies, V–V bond lengths and Mulliken spin densities for the ⁴A₂ and ⁶A₁ states of [V₄N₄]⁺. The energetic reference point is taken as the more stable of the two states.

Functional		4A2			6A1
	V–V/Å	ρ(V)	Erel/eV	V–V/Å	ρ(V)	Erel/eV
TPSS	2.64	0.86	+0.19	2.60	1.27	0
M06L	2.62	0.90	+0.38	2.56	1.29	0
BP86	2.65	0.85	0	2.61	1.28	+0.08
BP86-D3	2.65	0.85	0	2.61	1.28	+0.09
BLYP	2.68	0.85	0	2.63	1.28	+0.05
B3LYP	2.62	0.70	+0.65	2.62	1.36	0

), but the relative energies of the states are not: with the TPSS functional, the ⁶A₁ state is the global minimum, 0.19 eV lower in energy than ⁴A₂, while with the BP86 functional, the order is reversed, with ⁴A₂ 0.08 eV below ⁶A₁. These behaviours appear to be generally characteristic of the broad classes of functionals (Table 1): the meta-GGAs, TPSS and M06L [30] favour the ⁶A₁ state, while all three GGA functionals (BP86 [31,32], with and without dispersion corrections (Becke–Johnson D3 set [33]), and BLYP [34]) favour ⁴A₂. The hybrid B3LYP functional [35] strongly stabilizes the state of higher multiplicity, a trend that has been observed many times in the literature [36]. Clearly, the balance between one-electron terms, favouring double occupation of the most bonding orbitals, and electron–electron repulsions, favouring higher multiplicities, is a delicate one that challenges even those functionals that have previously proven suitable for metal–metal bonded systems (TPSS and BP86) [27].

In light of the conflicting conclusions about the relative stabilities of ⁴A₂ and ⁶A₁ that emerged from the DFT treatment, we have turned to multi-configurational approaches (specifically the Complete Active Space Self-Consistent Field (CASSCF) methodology) to offer an alternative perspective on the V–V bonding. In both cases, the calculations were performed using the geometry of the corresponding state obtained from the calculations with the TPSS functional, although the precise geometry of a given state is rather independent of the exact choice (with the exception of the M06L functional, which tends to afford rather shorter V–V bonds). The ⁴A₂ and ⁶A₁ states were described using a CAS(7,9) active space, where the nine orbitals include the 4e, 9a₁ and 12t₂ orbitals in Figure 2, and also the three additional orbitals of the 4t₁ set that are derived from the tangential $d_{xy}/d_{x^2-y^2}$ set. The natural orbitals (the orbitals that diagonalise the one-particle density matrix [37,38]) and their occupation numbers are shown in Figure 3 (the orbitals shown are for the ⁴A₂ state, but the iso-surface plots are almost indistinguishable from those for the ⁶A₁ alternative). The three t₂ orbitals had occupations close to 1.0, but the rather low occupations of the 4e orbitals (1.783) and the compensating high values for the 4t₁ set (0.143) are indicative of strong multi-configurational character. Indeed, the single configuration shown in the Kohn–Sham orbital diagram makes up only 78% of the ground-state wavefunction, the remainder being largely contributions from double excitations from 4e to 4t₁: these orbitals are bonding–antibonding combinations of the same set of basis functions, and so the strong correlations are consistent with Roos’s proposals for an appropriate choice of active space [39,40]. The natural orbital occupations for the ⁶A₁ state (shown in parentheses in Figure 3) confirm a rather similar picture of the bonding: the three components of the 12t₂ orbital were singly occupied, as were the 9a₁ orbital and one of the two components of 4e. The one remaining pair in the 4e orbital was strongly correlated with the three members of the 4t₁ orbital, but their overall occupations were reduced accordingly. In summary, the picture of bonding that emerged from the CASSCF calculations was qualitatively similar to the DFT perspective, with the caveat that the multi-configurational approach highlights the very strong static correlation within the doubly occupied V–V bonding orbitals. At the CASSCF level, the ⁴A₂ state was less stable than the ⁶A₁ alternative by 4.78 eV, but the incorporation of dynamic correlation using N-electron valence state perturbation (NEVPT2) theory reversed the order, such that the ⁴A₂ state was more stable by 5.88 eV.

2.2. Electronic Structure of [V₆N₆]⁺

Turning now to the larger [V₆N₆]⁺ cluster, we considered a range of structures that have been proposed in the literature for related species, including the columnar fragment of the rock-salt structure proposed by Hirabayashi and Ichihashi (D_2h symmetry) and also a hexagonal prismatic structure with two staggered V₃N₃ hexagons (D_3d symmetry), which has been identified as the global minimum for the Al₆N₆ cluster [41]. In the case of [V₆N₆]⁺, this hexagonal structure proved to be over 3 eV higher in energy than the D_2h alternative for all functionals, and was not considered further. As was the case for the [V₄N₄]⁺ cluster, we found a quartet ground state (⁴B_3g) with overall D_2h symmetry, with all alternative spin states substantially higher in energy. The V–V bond in the central layer of the cluster (V1–V1) is strikingly short at 2.23 Å (

Table 2. V–V bond lengths and Mulliken spin densities for the ⁴B_3g state of [V₆N₆]⁺.

Functional	4B3g
	V–V/Å 1	ρ(V) 2
TPSS	2.23, 2.52, 2.60	−0.10, 0.87
M06L	2.23, 2.51, 2.58	−0.12, 0.92
BP86	2.24, 2.53, 2.62	−0.03, 0.85
BLYP	2.27, 2.55, 2.64	−0.05, 0.86
B3LYP	2.24, 2.53, 2.59	−0.39, 1.07

¹ the three bond lengths correspond to V1–V1, V1–V2 and V2–V2 in Figure 4; ² spin densities are for V1 and V2.

), a value that is indicative of strong multiple bonding. The bonds between the central and upper/lower layers (V1–V2) were somewhat longer, at 2.52 Å, whilst those within the top and bottom layers (V2–V2) are longer still, at 2.60 Å. All V–V contacts are, however, shorter than those found in either the ⁴A₂ or ⁶A₁ states of [V₄N₄]⁺. The Kohn–Sham orbitals for the ⁴B_3g state, shown in Figure 4, confirm that, of the eleven available V 3d valence electrons, eight are paired in orbitals with significant V–V bonding character. Amongst the doubly occupied set, the 20a_g and 15b_3u orbitals are localized in the central V₂N₂ plane, and have strong V–V σ and V–V π character, respectively, consistent with the presence of a V=V double bond. The 4a_u and 10b_1g orbitals, in contrast, are bonding with respect to all eight V1–V2 contacts between the central plane and the upper/lower layers of the column, accounting for the relatively short distances of 2.51 Å. The remaining three valence electrons occupy the 21a_g, 14b_2u and 16b_3u orbitals, none of which are strongly V–V bonding or antibonding. There are, therefore, many qualitative parallels between the ⁴A₂ state of [V₄N₄]⁺ and the ⁴B_3g state of [V₆N₆]⁺: in both cases, some of the available electrons are paired up to form V–V covalent bonds, while the remainder are distributed singly over approximately non-bonding orbitals. The proportion of electrons in the two sets of orbitals changes as a function of the cluster size: for the ⁴A₂ state of [V₄N₄]⁺, four of the seven available electrons are involved in V–V bonding, while in the ⁴B_3g state of [V₆N₆]⁺ this fraction rises to 8/11, and the resulting increase in net bond order is behind the generally shorter bond lengths in the larger cluster.

The CASSCF wavefunction for the ⁴B_3g state of [V₆N₆]⁺ shown in Figure 5 revealed a qualitatively similar story about the V–V bonding, relative to the one established for [V₄N₄]⁺. We used a CAS(11,11) active space in which all four doubly occupied bonding orbitals in the Kohn–Sham scheme find antibonding counterparts in the virtual manifold that allow for strong correlations. These are most striking for the V–V σ and π bonds in the central V₂N₂ plane (occupations of the bonding–antibonding pairs are 1.726/0.275 (σ/σ*, 20a_g/13b_1u) and 1.705/0.317 (π/π*, 15b_3u/9b_2g)). Natural orbital occupations that differ dramatically from 2.0 and 0.0 are typical in clusters with multiple metal–metal bonds [42], and are highly indicative of weak bonding. The electrons in the 4a_u and 10b_1g orbitals that mediate the bonding between the layers are also strongly correlated (1.439/0.562 and 1.816/0.190 for the 10b_1g/11b_1g and 4a_u/5a_u pairs, respectively). The three remaining singly occupied orbitals are then 21a_g, 14b_2u and 16b_3u, all of which have occupations almost equal to 1.0. These three singly occupied orbitals bear a striking resemblance to the three singly occupied components of the t₂ orbital in the ⁴A₂ state of [V₄N₄]⁺.

3. Discussion

The two clusters considered here, [V₄N₄]⁺ and [V₆N₆]⁺, are both based on fragments of the rock-salt structure, as anticipated by Hirabayashi and Ichihashi in their original report of the mass-spectrometry experiments [14]. In the smaller cluster, quartet and sextet states are very close in energy, and meta-GGA functionals tend to favour the sextet, while GGA functionals predict the quartet to be the ground state. The B3LYP functional, as is typically the case, strongly stabilizes the higher spin state. Our analysis of the Kohn–Sham orbitals indicates that in the ⁴A₂ state, the cluster is stabilized by two pairs of electrons in the V–V bonding 4e orbital, giving a net V–V bond order of 1/3 and a perfectly tetrahedral geometry. The sextet (⁶A₁), in contrast, is D_2d-symmetric, with the distortion away from a perfect tetrahedron being driven by a first-order Jahn–Teller instability. The V–V bonds are shorter and stronger in the ⁶A₁ state as a result of the single occupation of the V–V bonding 9a₁ orbital, but this is compensated for by a weakening of the V–N bonds. CASSCF calculations (CAS(9,7) active space) are consistent with the picture of bonding that emerges from the DFT calculations, but also reveal the strongly multi-configurational nature of the V–V bonding, which leads to low (<<2.0) occupancies of the 4e orbitals with compensating populations in the V–V antibonding 4t₁ set. The ground state of the larger [V₆N₆]⁺ cluster appears to be less controversial: a ⁴B_3g state is the most stable for all functionals considered here, and also at the NEVPT2 level. The cluster is stabilized by a strong V=V double bond between the metal atoms in the central V₂N₂ layer, a feature that emerges irrespective of the chosen theoretical model. In the DFT calculations, the 20a_g (V–V σ) and 15b_3u (V–V π) orbitals are doubly occupied, and the three unpaired electrons in the ⁴B_3g ground state are localized on the upper and lower layers. To a first approximation, the cluster can therefore be viewed as a diamagnetic central V₂N₂ unit, with radical character delocalized over the upper and lower planes. CASSCF calculations support this basic model, but the multi-configurational character of the V=V double bond is again prominent, with occupations of the σ and π bonding and antibonding orbitals around 1.70 and 0.30. The balance between localized multiple V–V bonding in [V₆N₆]⁺ and the more delocalized situation in [V₄N₄]⁺ hints at a subtle relationship between structural and electronic properties, which may have important consequences for the reactivity of the nitrogen-deficient clusters, [V_n+1N_n]⁺. These will be the subject of subsequent studies.

4. Computational Methods

All DFT calculations reported in this paper were performed using Gaussian 16, revision A.03 [43] while the CASSCF calculations were done with ORCA, version 4.0 [44]. Density functional calculations were performed using Ahlrichs’s def2-TZVP basis set [45] on V and N either using the meta-GGA TPSS or M06L functionals, the GGA alternatives, BLYP or BP86, or the hybrid B3LYP. CASSCF calculations were performed using the same TZVP basis set and the RIJCOSX approximation to the exchange terms [46]. Dynamic correlation was accounted for using the N-electron valence state perturbation theory (NEVPT2) [47]. All optimised energies and geometries are summarised in Table S1, Supplementary Materials.

5. Conclusions

In this paper we used a combination of density functional theory and correlated ab initio theory (CASSCF) to explore the nature of the bonding in two clusters, [V₄N₄]⁺ and [V₆N₆]⁺, which have been observed in mass-spectrometry experiments. Unlike clusters with similar composition, both [V₄N₄]⁺ and [V₆N₆]⁺ are remarkably resistant to further reaction with NH₃, suggesting that they share some unique structural and/or electronic features. Our calculations confirm that the two clusters are based on cubic architectures and can be viewed, to a first approximation, as fragments of the rock-salt lattice. However, it is clear that direct covalent V–V bonding plays a very substantial part in stabilizing these structures, just as it does in the low-valent solid phases with similar compositions, such as VN and VO. This bonding is particularly well-developed in the quartet ground state of the [V₆N₆]⁺ cluster, where a localized V=V double bond is present in the waist of the cluster. The DFT and CASSCF approaches offer broadly equivalent pictures of the metal–metal bonding, but the multi-configurational approach highlights the very strong static correlation that is typical of metal–metal bonds, particularly amongst the first transition series. The strong V–V bonds, along with the very symmetric structure, are consistent with the observed resistance of these two clusters to further reactions. Future work will focus on the interplay between V–V bonding and catalytic activity for the less symmetric nitrogen-deficient clusters, [V₄N₃]⁺ and [V₆N₅]⁺, both of which are substantially more reactive [14].