Charged vanadium-benzene multidecker clusters: DFT and quantum Monte Carlo study

One-dimensional transition metal-benzene sandwich clusters are among the most important organometallics^1–3 with a range of applications in catalysis, polymers, molecular recognition,⁴ and high-density storage. However, perhaps the most intriguing applications of such one-dimensional nanowire systems are in spintronics⁵ and quantum computing,⁶ where different spin orientations have unique transport properties. Indeed, several of transition metal organometallic sandwich structures have been predicted to be half-metallic ferromagnets,^7–14 predominately by density functional theory (DFT) studies. However, the more accurate quantum Monte Carlo (QMC) calculations¹⁵ found these systems to be ferromagnetic insulators with large and broadly similar spin gaps. Another important area of research is the behavior of magnetic adatoms on graphene; benzene-transition metal half-sandwiches have often been used as models for such nanostructures.¹⁶ QMC methods provided new insights that have been markedly different from DFT,^17,18 indicating the importance of a proper treatment of electronic correlation in these systems.

Experimentally, the information on these systems is often gained from studies of their charged counterparts, cations, and anions.¹⁹ However, interpretation of such experiments with regard to neutral systems might be difficult or even lead to biased results.¹⁷ This makes accurate theoretical predictions on both neutral and charged complexes particularly valuable since it would be directly comparable to experiments as well as shed new light on the differences between the neutral and charged cases.

As discussed above, the electronic correlation of transition metal nanostructures makes these systems difficult to describe at a single-particle level even in neutral state. An even more complicated situation arises with anions, where the majority of commonly used DFT exchange-correlation (xc) potentials describe anions as having only a fraction of the extra electron bound, while the remaining fraction drifts off to infinity when sufficiently flexible basis sets are employed. The common strategy to alleviate this deficiency is to employ long-range corrected (LC) xc functionals along with flexible basis sets with diffuse functions.²⁰ It is not uncommon that such an approach can lead to unbound anions where experiments strongly suggest stable anions. On the other hand, the use of generalized gradient approximation (GGA) functionals yields stable anions in very good agreement with experimental electron affinities (EAs). Unfortunately, detailed analysis shows that this favorable agreement depends very significantly on error cancellations. We find that truly correct stable anions result only if LC xc functionals are used for generating the initial orbitals that enter into explicitly correlated wave function methods, such as QMC. Anions are also challenging in many-body wave function methods since electron correlation has to be recovered at a much higher level of accuracy than for most neutral systems. In fact, only a handful of QMC calculations on anions have been performed.^21,22 On the contrary, the more strongly bound cations are comparatively easily described by the mainstream GGA-DFT methods.

We demonstrate these findings on some of the most studied transition metal-benzene systems, namely, charged vanadium-benzene multidecker clusters V_nBz_n+1 for n = 1–3. These structures can be seen in Fig. 1. We determine the stable spin multiplets, fragmentation energies, ionization potentials, and electron affinities. In addition, we will also briefly consider the charged half-sandwiches (VBz).

The presented calculations used the following four-level refinement strategy which enabled us to eliminate basically all systematic biases: (1) initial geometries were obtained from DFT optimization, (2) trial wave functions were constructed from a DFT calculation, (3) trial wave functions were optimized using VMC (variational Monte Carlo) techniques, and (4) energies were computed from DMC (diffusion Monte Carlo) simulations. For DFT modeling, we used the GAMESS suite of codes,²³ while all VMC and DMC calculations used the QWalk code.²⁴ 

The ground-state geometries were initially obtained using DFT techniques spanning a number of xc-functionals. We have used the Becke-Perdew-Wang 91 (BPW91)^25,26 as a representative GGA xc-functional, Becke-three-parameter-Lee-Yang-Parr (B3LYP)²⁷ as a representative hybrid xc-functional, and Coulomb-Attenuating Method B3LYP (CAMB3LYP) functional²⁸ as a representative LC xc-functional. The same xc-functionals were also used in the DFT energy calculations and in construction of trial nodal hypersurfaces using determinants with DFT orbitals. In addition to the systematic construction of trial wave functions, use of the BPW91 functional was motivated by the fact that previous DFT studies of transition metal-benzene systems all used the BPW91 functional^29–31 or the related Perdew-Burke-Ernzerhof (PBE)^32,33 functional^9,16 that enables us to make direct comparisons with previous DFT calculations.^9,29–31 The Greeff-Lester type of effective core pseudopotential^34,35 and cc-pVTZ and aug-cc-pVTZ³⁶ basis sets were used for all species. Casula T-moves algorithm for variational evaluation of non-local pseudopotentials³⁷ was used. We note that we found the DFT calculations with diffuse aug-cc-pVTZ basis sets very hard to converge. For that reason, a complete systematic study of all multiplets with the aug-cc-pVTZ basis was only performed for the V₂Bz₃ anion; for all other systems, only the respective ground-state multiplets were studied. Dynamical correlations were explicitly built in via Schmidt-Moskowitz Jastrow factor,^38,39 including electron-electron, electron-nucleus, and electron-electron-nucleus correlations. A single determinant of DFT states was used to determine the nodal part of the trial wave function with parameters of the Jastrow factor optimized. Hence, keeping in mind the limitations of the previous DFT studies, accurate QMC calculations for these systems have been performed in parallel to the DFT modeling.

We have carefully checked the effect of structural relaxation due to the above mentioned subtleties, especially on the anions. As will be shown below, the most difficult anion to describe is the $V2Bz3−$⁠. For this system, we compared the DFT geometries and found only negligible structural differences. In addition, when recomputed with a different xc-functional, the resulting energy differences showed only very minor variations. We attempted a QMC structural optimization focusing only on the transition metal atom-benzene distance.^15,17 Again, the atomic structure-related differences were significantly smaller than the chemical accuracy. Hence, all the results were obtained using the DFT structural optimizations.

In order to correctly describe bonding of the species with an electron added or removed, we first study the electron affinities and ionization potentials. We start with the anions and use the most difficult system, V₂Bz₃, as an example described below. Only the respective ground-state multiplets, see Sec. III B, will be considered.

The results for both adiabatic EAs and vertical detachment energies (VDEs) for all the functionals (BPW91, B3LYP, and CAMB3LYP) and basis sets studied (cc-pVTZ and aug-cc-pVTZ) are shown in Table I; all calculated DFT/QMC energies used to determine the energy differences can be found in the supplementary material.⁴⁰ Surprisingly, the simple GGA BPW91 xc-functional yields remarkable agreement with experiment. This finding is consistent with the work of Masabuchi et al.,¹³ who found similar results with a GGA functional. Such agreement with experiment is surprising since the localized d electrons are unlikely described correctly with an ordinary GGA functional. We note that if more accurate QMC electronic correlations are built using the trial wave functions with BPW91 orbitals, both EA and VDE yield unstable anions, see Table I. This suggests that the agreement with experiment is indeed due to fortuitous error cancellations. As shown in  Appendix B, the most obvious failure of the BPW91 functional is its inability to correctly describe the spin structure of the anions. We have tested a series of improvements to GGA at the DFT level: (1) localizing d electrons by use of the B3LYP hybrid functional; (2) B3LYP with basis sets improved by diffuse bases (denoted B3LYP^⋆ in Table I and henceforth); (3) the LC xc-functional CAMB3LYP; and (4) CAMB3LYP with basis sets improved by diffuse bases (CAMB3LYP^⋆ in Table I). From Table I we conclude that regardless of the treatment, the EA/VDE remains negative/barely positive; this renders the anions unstable/barely stable, which is in disagreement with the experimental electron affinity of 0.35 ± 0.05 eV. In order to illustrate the trends, the VDE DFT results are compared to experiment also in Fig. 2. Hence, for our system, use of LC DFT functionals along with diffuse bases is insufficient to correctly describe the anions, and a more accurate treatment of electronic correlation is required.

Clearly, these DFT results leave some questions unanswered and the problem calls for more accurate approaches. QMC calculations of anions have been carried out on a handful of systems in the past; however, most of the applications were devoted to atomic or rather small molecular systems.^42,43 Since the binding of an additional electron to a neutral system is predominantly a correlation effect, QMC methods can provide both accurate predictions and also valuable insights. The QMC applications to systems with first- or second-row elements were quite successful since the relevant states were either s- or p-based in the outermost shell; the main challenge was describing the subtle medium- and long-range correlations. Here, QMC is very effective and produces reliable estimations of related energetics since nodal errors are secondary and medium-/long-range correlations are described very accurately.

For 3d-elements the additional electron usually goes into the open 3d-subshell; however, that is not always the case. There are exceptions such as Sc⁻ where the additional electron goes into the 4p-subshell and the resulting state is a singlet, i.e., the state is an exception to the Hund’s rule.⁴³ The correlation effects are very significant and their treatment require very high level of many-body treatment. For example, calculations of the electron affinity of the Fe atom (neutral ⁵D(3d⁶4s²)) are challenging since the additional electron goes into the d-shell (anion ⁴F(3d⁷4s²)) that is rather strongly localized.⁴² The anionic state mixes a large number of excitations with sizable exchange and correlation contributions; however, these are counteracted by the strong electron-electron (on-site Hubbard U) repulsion so that the resulting affinity is very small ≈0.15 eV.⁴² Note that the competition between Coulomb repulsion and exchange with correlation can easily tilt the result to a nonbinding (negative) affinity, e.g., the Mn atom. Accurate calculations of such small energy differences with large accompanying correlation energies (in the 3s3p3d4s space of the order of ≈20 eV) require highly accurate trial wave functions with a high quality correlation description.⁴³ Capturing all of the relevant competing effects is, in general, very demanding.

To some extent, these effects are also present in our systems but are further complicated by the fact that the anions involve both short-range correlations in the vanadium d-shell as well as long-range correlations in the tails of anionic states. The calculations show that the accuracy of the trial function plays an important role; the balance between short-range and long-range correlations has to be captured in a systematic manner.

The picture that has emerged from using QMC methods for these effects is shown in Table I and Fig. 2. While both BPW91 and B3LYP nodal hypersurfaces yield negative EA and VDE values, adding diffuse basis sets to B3LYP accounts for a qualitative improvement of the QMC results. Use of the LC functional (CAMB3LYP) and diffuse basis sets further improves the results with the VDE error bars overlapping with the experimental error bars. The effect of diffuse basis sets appears larger than the improvement from a more sophisticated xc-potential. Unlike ionization potentials, where experiments clearly suggest that a correspondence to vertical ionization potentials (VIPs),⁴¹ the experiments for electron affinities¹³ do not specify a correspondence to EAs or VDEs; however, the accompanying GGA-DFT calculations quote only EAs. Obviously, our results correspond more closely with VDEs, as seen in Fig. 2.

From the analysis above, we conclude that in our calculations, only LC DFT functionals with appropriate diffuse basis sets included (CAMB3LYP^⋆) provide nodal hypersurfaces accurate enough to stabilize the V_nBz_n+1 anions. Hence, all anionic calculations have been performed in CAMB3LYP^⋆ approximation. The EA/VDE results for the entire V_nBz_n+1 n = 1–3 series along with the results for the VBz half-sandwich structure are shown in Table II using DFT/BPW91, DFT/CAMB3LYP^⋆, QMC/BPW91, and QMC/CAMB3LYP^⋆ approximations. For the sake of showing the trends, our VDE results are also shown in Fig. 3. On average, VDEs are again in better agreement with the experiments than the EAs. The DFT/BPW91 results for the V_nBz_n+1 n = 1–3 agree nicely with the measured electron affinities. The DFT/CAMB3LYP^⋆ is a clear improvement over the GGA results. The stability of the V_nBz_n+1 n = 1–3 series is lower than in the experiment, with V₂Bz₃ being strongly affected. The impact from the inclusion of more accurate electronic correlation in the QMC/CAMB3LYP^⋆ treatment is clearly visible; the VDEs of the $V3Bz4−$ are within chemical accuracy with respect to the experiment and agree qualitatively for $VBz2−$⁠, where only the negative sign could be determined experimentally. The deviation from the experimental value for $V2Bz3−$ is marginally larger. For the VBz⁻ half-sandwich, the behavior is qualitatively similar to that of the full sandwiches. The DFT/BPW91 approach yields excellent agreement with experiments. However, the inclusion of correlation at the QMC/BPW91 level yields a VDE that is ≈0.25 eV too high, highlighting the fact that the BPW91 nodal hypersurfaces are inadequate for these effects as well as for the description of the xc-hole. Only modeling the VBz⁻ at the CAMB3LYP^⋆ level brings both DFT and QMC results into the agreement with the experimental value.

The calculated results for VIPs for the V_nBz_n+1 n = 1–3 series and for the VBz half-sandwich are compared with the experimental results in Table III and in Fig. 4. The more tightly bound cations are very well described by the GGA BPW91 functional for all systems considered, and the QMC correlations produce only a minor improvement for the full sandwiches and no improvement for the half-sandwich. For that reason, DFT/BPW91 and QMC/BPW91 approaches will be used in all other calculations for cations.

The calculated relative energies of the $VnBzn+1−$ anions for n = 1–3 are shown in Fig. 5. Note that all of the CAMB3LYP^⋆ spin multiplets were calculated for $V2Bz3−$ and only the minimum energy multiplets for $VBz2−$ and $V3Bz4−$⁠. All the energy differences are in excess of ≈0.5 eV; these differences are well reproduced at the DFT/BPW91 and QMC/BPW91 levels. The ground-state $VnBzn+1−$ multiplets, n = 1–3, are equal to n, the number of vanadium atoms in the cluster. Compared to the neutrals, which had the vanadium spins ferromagnetically coupled and followed the n + 1 rule,¹⁵ the multiplicity of the anion is decreased by one, i.e., have one spin quenched. This rule can clearly be seen in the calculated spin densities, shown in Fig. 5. The DFT and QMC results for the half-sandwich anion, shown in Fig. 5, indicate that DFT and QMC find different ground-state multiplets, namely, a quintet in DFT/CAMB3LYP^⋆ and a triplet in QMC/CAMB3LYP^⋆.

The calculated relative energies of the $VnBzn+1+$ cations for n = 1–3 are shown in Fig. 6. The ground-state $VnBzn+1+$ multiplets, n = 1–3, are equal to n + 2, the number of vanadium atoms in the cluster plus 2. Compared to the neutrals, which had the vanadium spins ferromagnetically coupled and followed the n + 1 rule,¹⁵ the vanadium atoms are also ferromagnetically coupled, as can be seen from the calculated spin densities shown in Fig. 6. The DFT and QMC results for the half-sandwich cation in Fig. 6 indicate that also for the cation, DFT and QMC find different ground-state multiplets, namely, a triplet in DFT and a quintet in QMC.

Once we identified the correct ground states and their spin multiplets, we determined the energetically preferred fragmentation of the charged vanadium-benzene species into vanadium atoms and benzene molecules and their respective ions. The results for the anions are shown in Table V. In the following, we label the multiplicity, M, or atomic/molecular state in brackets. In all cases, the lowest energy fragmentation channel is the dissociation into neutrals and the V⁻ anion,

which requires ≈3.4-3.6 eV/V atom. Fragmentation into neutral species,

is only very marginally less favorable (≈0.15 eV (n = 3) to 0.4 eV (n = 1)) than the dissociation into V⁻. This energy ordering of fragmentation channels is interesting, given the very small values of VDE, especially for V₂Bz₃. The lowest energy fragmentation into charged benzene species,

requires between ≈0.6 eV (n = 3) and ≈1.8 eV (n = 1) more energy. Energetics of the half-sandwich anion fragmentation follows the same principle, as shown in Table IV. The impact of electronic correlation is clearly seen in the DFT fragmentation energies where the CAMB3LYP* consistently underbinds by about 1 eV with respect to the QMC results; for the BPW91 results, see Table SXI of the supplementary material.⁴⁰ 

The results for cations are summarized in Table V. Fragmentation into vanadium cation,

is strongly energetically favored on both DFT and QMC levels by ≈1-2.5 eV over the fragmentation into charged benzene molecules,

Again, the energetics of the half-sandwich cation fragmentation follows the same principle, as shown in Table V. DFT BPW91 overbinds by about ≈0.4 eV with respect to the QMC results.

We have presented a study of charged vanadium-benzene half- and full-sandwich systems, including both anions and cations. The study is based on the use of explicitly correlated QMC methods with the only essential approximation being the fixed-node approximation. These highly correlated calculations allow for direct comparisons with the various DFT approaches which we use to fix the nodal hypersurfaces. In particular, we have studied the ground-state spin multiplets of $VnBzn+1±$ in the range of n = 1–3 and of the VBz^± half-sandwiches. We have determined their ionization potentials and electron affinities, as well as the lowest-energy fragmentation channels.

We show that the electronic correlation in these systems plays a crucial role, especially in the anions. Commonly used DFT xc-functionals, e.g., gradient corrected, hybrid, and range-separated hybrid functionals, fail to accurately describe these systems. Being a single determinant theory, it is frequently found that the low-spin DFT calculations on these systems show a large amount of spin contamination, a feature that would be, via the nodal hypersurfaces, transferred also to the DMC calculations. Fortunately, only mild spin contamination is found in our DFT wave functions, see Tables SIII-SVII of the supplementary material.⁴⁰ The most critical system we study is the V₂Bz₃ anion, where only a combination of the LC-corrected DFT theory for nodal hypersurfaces with electronic correlation built in through QMC methods yields stable anions in agreement with the experimental findings. The cations are much more easily modeled by DFT techniques, where QMC methods provide mainly quantitative improvements. The exceptions are the ground-state multiplets of the VBz^± half-sandwiches, where the correlated QMC treatment provides different ground-states than the DFT methods and the well-known overbinding (GGA)/underbinding (hybrids) in calculation of the dissociation energies. Surprisingly, the influence of the atomic structure on the resulting energies was found negligible in the different treatments with hybrid functionals, and the differences in computed energies result almost entirely from electronic correlation effects. On the other hand, the different spin structures in GGA and hybrid functionals, Fig. 8, lead to more significant differences also in their atomic structures, for additional examples see Fig. S2 of the supplementary material.⁴⁰ The size of the systems we have been able to treat with the correlated QMC methods significantly raises the standards and expectations of quantitative modeling of electronic structure using stochastic methods.

This research was supported in part by No. APVV-0207-11 and VEGA (Nos. 2/0007/12 and 2/0162/15) projects. L.M. gratefully acknowledges support by the NSF Grant No. DMR-1410639. We also gratefully acknowledge use of the Hitachi SR16000/M1 supercomputer system rat CCMS/IMR, Tohoku University, Japan and XSEDE allocation at TACC.

In Sec. III A, we have argued that the different DFT xc-functionals found the anions to be unstable/barely stable depending on the employed xc-functional. Only LC xc-functionals used along with diffuse bases gave qualitatively correct DFT results across the entire series of studied systems; quantitatively correct results could only be obtained after the inclusion of QMC electronic correlations. The most striking example is the $V2Bz3−$ anionic sandwich structure which most of the functionals, contrary to the experiments, modeled as unstable. In order to better understand those differences, we directly compared the different underlying nodal hypersurfaces (BPW91, B3LYP, B3LYP^⋆, CAMB3LYP, and CAMB3LYP^⋆). We note that the QMC energies are critically dependent on the nodal hypersurfaces.³⁹ To this end, we use the procedure of direct stochastic “measuring the differences between nodal hypersurfaces,” ΔS_node, as outlined in Ref. 49. Similarity between a pair of nodal surfaces can be quantified using the value of coefficient of the fall-off ξ₁ in the fitting function, $y=ξ0exp(−ξ1x)/x3$⁠,⁴⁹ used to approximate the difference of the distributions of the nodal hypersurfaces. The distribution of differences between the different DFT nodal hypersurfaces, ΔS_node, is shown in Fig. 7. We indeed confirm the trends of fairly tiny differences between BPW91 and B3LYP (large ξ₁), and larger differences resulting from adding the diffuse bases and use of LC-xc functionals (smaller values of ξ₁).

The most striking difference between the different DFT approaches for the description of anions was the difference between the BPW91 and the hybrid functionals, shown in Fig. 2. The DFT BPW91 description resulted in VDEs which were essentially in perfect agreement with the experiments, while all the hybrid approaches resulted in an unstable anion. However, when those orbitals were used to construct the nodal hypersurfaces for the QMC calculations, this trend was reversed and the hybrid functionals resulted in a more correct description of the multidecker anionic clusters.

We show here that the main difference comes from significant deviations in spin densities resulting from the two types of xc-hole models. In fact, this difference is also found in the atomic structure of the $V3Bz4−$ anion, which exhibits a bended geometry in the BPW91 but a straight structure in the CAMB3LYP^⋆ description shown in Fig. 8. The spin densities obtained in the other DFT treatments (B3LYP, B3LYP^⋆, and CAMB3LYP) are qualitatively similar to the CAMB3LYP^⋆, as shown in Fig. S2 of the supplementary material.⁴⁰ As can be seen from Table VI, while the spin contaminations in BPW91 and CAMB3LYP^⋆ are very similar, the resulting spin structures are completely different. The CAMB3LYP^⋆ generates a spin density corresponding to quenching of one spin of the ferromagnetically coupled spins on the V atoms upon adding an electron, in complete agreement with the expected behavior. However, the BPW91 xc-functional features a quenched spin density only for $VBz2−$⁠, while for $V2Bz3−$ and $V3Bz4−$⁠, no vanadium spin is completely quenched. This is due to the over-delocalized nature of the vanadium d electrons in the BPW91 model, in which the added electron avoids occupation of a localized atomic-like 3d state. On the contrary, localization of the added electron into the atomic-like 3d state is precisely the effect the hybrid functionals are able to achieve. We presume that this feature is responsible for the fact that the DFT VDE values are apparently correct while the electronic structure and the underlying nodal hypersurfaces are incorrect.