First principles studies of the electronic and structural properties of the rutile VO2(110) surface and its oxygen-rich terminations

We present a density functional theory (DFT) study of the structural and electronic properties of the bare metallic rutile VO2 (110) surface and its oxygen-rich terminations. Due to the polyvalent nature of vanadium and abundance of oxide phases, the modelling of this material on the DFT level remains a challenging task. We discuss the performance of various DFT functionals, including PBE, PBE + U (U = 2 eV), SCAN and SCAN + rVV functionals with non-magnetic and ferromagnetic spin ordering, and show that the calculated phase stabilities depend on the chosen functional. We predict the presence of a ring-like termination that is electronically and structurally related to an insulating V2O5 (001) monolayer and shows a higher stability than pure oxygen adsorption phases. Our results show that employing the spin-polarized SCAN functional offers a good compromise, as it offers both a reasonable description of the structural and electronic properties of the rutile VO2 bulk phase and the enthalpy of formation for oxygen rich vanadium phases present at the surface.


Introduction
Vanadium dioxide undergoes a metal-to-insulator transition (MIT) at 68 °C which makes it usable in many applications such as optoelectronic switches, sensors or Mott-field effect transistors [1,2,3].Below the transition point, VO 2 undergoes the phase transition from a metallic rutile (R) phase (space group P4 2 /mnm) to a semiconducting monoclinic (M) phase (space group P2 1 /c), forming chains of paired vanadium atoms and thereby diminishing its electrical conductivity.The rutile phase is paramagnetic above the transition temperature while the monoclinic phase is considered to be nonmagnetic [4].The transition temperature can be driven towards room temperature by doping with high valence cations such as W (6+) [5] to make this material available in real applications, such as smart window coatings.
Owing to the numerous potential applications of VO 2 , the bulk properties have been studied applying a vast number of theoretical approaches [6].It has been shown that standard methods within the DFT framework fail to completely recover the electronic, structural and magnetic properties of the rutile and monoclinic VO 2 phase.For example, both the local density approximation (LDA) and the generalized gradient approximation (GGA) fail to open an 0.6 eV band gap in the monoclinic (M1) phase [7,8].This shortcoming might be eliminated either by using the DFT+U variant [9], or meta-GGA and hybrid functionals, as analyzed by Stahl and Bredow [10,11].According to their findings, DFT+U leads to strong structural distortions of the monoclinic phase and a wrong energetic ordering of the VO 2 phases.Furthermore, conventional hybrid functionals such as HSE06 or PBE0 lead to a splitting of the conduction band in the metallic rutile phase and thus the Fock mixing parameter needs to be adjusted to correctly describe the electronic structure of the rutile and monoclinic phases [12,11].Stahl and Bredow also concluded [10] that the meta-GGA SCAN functional offers a good compromise between accuracy and computational cost.
The driving mechanism behind the MIT was recently investigated by Brito et al. [13] by performing combined DFT and embedded dynamical mean-field theory (DMFT) calculations.Such a combination constitutes a powerful approach to evaluate the many-body electronic structure of strongly correlated materials [14] and in particular to describe the MIT between a metal and a Mott insulator, as is the case for the present system, where the electronic transition can be described as a Mott transition in the presence of strong intersite exchange.This is in agreement with the finding of Zhu et al. [6] who used the modified Becke-Johnson exchange together with LDA correlation potentials for both rutile and monoclinic VO 2 phases, to show that the MIT can be characterized as a correlation driven transition.
In contrast to VO 2 bulk systems, VO 2 surfaces and VO 2 thin films have been studied to a much lesser degree.The surface free energy of bare and oxygencovered low-index facets was investigated by Mellan et al. [15] employing the non-magnetic PBE functional.PBE+U calculations performed by Wahila et al. [16], indicate that the surface free energy of corresponding surfaces is lower for the rutile than for the monoclinic phase.Furthermore, they find it necessary to include spin-polarization in order to avoid negative values for the surface free energies of oxygen-rich reconstructions.In a very recent study [17], the rutile and monoclinic VO 2 surfaces were studied with the hybrid sc-PBE0 functional with an adjusted Fock mixing parameter.Despite its good performance on both bulk phases, the hybrid sc-PBE0 functional surprisingly fails to describe the surface energies of rutile VO 2 which do not converge upon increasing the thickness of the slab and even yield negative values for thicker slabs.On the experimental side, recent studies show the presence of a (2 × 2) reconstruction on the VO 2 (110) surface, accompanied by a change in the surface oxygen content of single crystals [18] and thin films [19,16].
The aforementioned DFT studies using the PBE(+U) approach on the rutile VO 2 (110) surfaces under oxygen-rich conditions [16,15] were done only for regular oxygen adsorption phases.However, a very recent study by Wagner et al. [18] indicates the presence of (2 × 2) surface reconstructions with tetrahedrally coordinated V atoms.In order to shed some more light on these findings, we present a detailed study of various oxygen adsorption phases and surface reconstructions on rutile VO 2 surfaces, and evaluate the performance of common DFT functionals, including GGA, GGA+U,meta-GGA and meta-GGA+U.In agreement with the experimental findings of Wagner et al. [18], we show that the simple regular oxygen adsorption phases considered in the previous calculations [15,16] are less stable than surface terminations with tetrahedral V coordination polyhedra that are structurally and electronically related to a V 2 O 5 (001) monolayer.This findings impose another important constraint to the DFT functional used for the description of off-stoichiometric VO 2 surface terminations, namely a correct description of the stability of the insulating V 2 O 5 phase with respect to VO 2 .

Computational details
All calculations were performed with the Vienna Abinitio Simulation Package (VASP) [20,21].We used the projector augmented wave (PAW) method [22,23] for treating the core electrons.The Bloch-functions for the 6 valence electrons of oxygen (2s 2 2p 4 ) and the 13 for vanadium (3s 2 3p 6 3d 4 4s 1 ), were expanded in a plane wave basis set with an energy cut-off 500 eV for all final relaxations.The Brillouin zone was sampled with the Γ-centered Monkhorst-Pack scheme [24] using 1.4 k-points/ Å−1 to ensure the absolute convergence within the accuracy 1 meV per atom.The electronic optimization self-consistent loop was converged to the 10 −5 eV and ionic relaxation was stopped when residual forces were smaller than 10 −2 eV A −1 .
The local electronic structure was also studied with help of the Wannier90 code [25].All sub-bands that comprise O 2p and V 3d bands were projected onto the same amount of Wannier functions.The local coordinate system of the Wannier functions localized at vanadium atoms was rotated by 45°around the z−axis as proposed by Eyert [7] in all systems studied.This rotation principally also implies an exchange of d xy and d x 2 −y 2 orbitals, but in our work we keep the naming of the orbitals the same as in the original basis.
We used several DFT functionals to describe the exchange-correlation energy, namely the PBE, PBE+U, SCAN [26], SCAN+U and SCAN+rVV [27].The influence of the Hubbard U within Dudarev's implementation [28] on the properties of VO 2 phases has been studied in the work of Stahl et al [10], showing that values of U between 1.4 eV to 2.6 eV correctly open a band gap in the monoclinic phase, but not in the rutile phase if the experimental structures are used in the calculations.We therefore settle on an U value of 2 eV as used before in point defect calculations [29].Based on these functionals we performed both nonspin polarized (NM) and spin-polarized calculations of ferromagnetic arrangements (FM) which turned out to be energetically most favourable for the rutile phase [10].
All slab calculations were performed at bulk optimized lattice constants applying the corresponding DFT functional and spin ordering.The stability of the respective VO 2 terminations was determined via evaluation of the surface free energy, calculated by linear regression from total energies of 5 to 8-layered slabs, as described by the following scheme.The surface energy of a slab that is composed of N bulk formula units, having a surface area S is extracted from DFT calculations according to the known formula: where E slab and E bulk are the total slab and bulk energies obtained from DFT calculations.A direct evaluation of σ requires a different sampling of Brillouin zones for bulk and slab models.In order to avoid this drawback, this formula can be rearranged towards a slab energy as a function of the number of layers, Comparison of these equations yields an expression for the surface energy as: where q is a parameter obtained from a linear regression based on calculations for several slab thicknesses.
To calculate the surface free energy of the most stable off-stoichiometric terminations we follow the procedure described in the work of Reuter [30].We approximated the Gibbs free energy with the total DFT energies of the slabs and reference systems in the oxygen-rich limit so that all vibrational contributions are neglected.The formula used to extract the surface free energy from a slab that consists of N V vanadium atoms and N O oxygen atoms as a function of oxygen chemical potential µ O reads where E VO2 and E O2 are reference energies for the total energy of the rutile VO 2 bulk and oxygen molecule respectively, E slab marks the total energy of the slab with the surface area S. We used 5-layered slabs with the symmetric top and bottom surfaces with a separating vacuum layer kept at 15 Å to study all (2 × 2) surface reconstructions.The values of the oxygen chemical potential are considered for a range of energies where the VO 2 phase is thermodynamically stable.Its oxygen-rich limit is approximately given by the reaction enthalpy of the oxidation reaction of the VO 2 phase: The experimental value −1.28 eV is determined from the formation enthalpies of the VO 2 [31] and V 2 O 5 [32] while the calculated value is considered in the 0 K temperature limit by neglecting vibrational and entropic contributions.This approximation yields an 1.3 % error in experimental heat of formation of vanadium pentoxide compared to the value obtained at the room temperature [32].The reaction enthalpyis therefore calculated using the total energies of bulk systems and oxygen molecule as In order to sample the large configuration space of the VO 2 (110) (2 × 2) terminations, we performed an optimization of the random structures generated by the USPEX package [33,34,35].A more detailed information on this optimization procedure can be found in the supplementary material.

Bulk properties of the VO 2 (R) and V 2 O 5 phases
The high-temperature rutile VO 2 phase exhibits a tetragonal structure (space group P4 2 /mnm, a=4.55 Å, c=2.86 Å) [36] with two VO 2 units that form the unit cell.Oxygen atoms are arranged in octahedra around vanadium atoms with a slightly shortened z-axis that is arranged perpendicular to the c-lattice vector as visualized in Figure 1 (right panel).In contrast, the V 2 O 5 phase shows an orthorhombic structure (space group P mmn , a=11.51 Å, b=3.56 Å, c=4.37 Å) [37] that is composed of van der Waals bonded layers stacked along the [001] direction.Note that there is certain ambiguity in the notation of the lattice vectors and in some papers the b-and c-lattice vectors are interchanged.In the present work we consider the V 2 O 5 layers to be oriented parallel to the (001) plane as shown in the left panel of Figure 1.
While the structural and electronic properties of the rutile VO 2 phase were calculated using both ferromagnetic and non-magnetic spin configurations, only closed-shell calculations were performed for the V 2 O 5 phase because no unpaired electrons are present.
Both phases were treated with all the aforementioned DFT functionals.The results are summarized in Table 1.Considering the rutile VO 2 phase, the lattice parameters are found to be well described well both with spin-polarized PBE and SCAN(+rVV) functionals, leading to errors below 1 %, while PBE+U and SCAN+U functionals tend to overestimate the c-lattice constant by 6 %, causing a similar overestimation of the c/a ratio.On the other side, non-magnetic calculations underestimate the clattice vector approximately by ∼3 % while the alattice vector turns out to be either overestimated by 2 % (PBE and PBE+U) or in good agreement with the experimental values (SCAN and SCAN+rVV).Consequently, the c/a ratio is calculated as 3-5 % lower than the experimental value.
Concerning the V 2 O 5 phase, all functionals yield in-plane a-and b-lattice parameters in very good agreement with experiment, with relative errors below 1.5 %.
For the c-lattice vector perpendicular to the layers the agreement is worse, being either overestimated by 11 % for PBE and PBE+U, or underestimated by 3 % for SCAN, by 4 % for SCAN+U and 5 % for SCAN+rVV.An overestimation by 6 % to 8 % in case of the PBE functional has already been reported in previous studies [38,39,40].According to our calculations, a slightly enlarged interlayer distance is accompanied only by a small energy penalty and thus our values differ slightly from the reported values while the a-and b-lattice constants match very well.It should also be noted that the spacing between the V 2 O 5 layers is already underestimated on the SCAN level, and that the additional van der Waals contributions in the SCAN+rVV functional lead to an even smaller value.Concerning the unit cell volumes of rutile VO 2 , non spin-polarized calculations as expected deliver slightly smaller volumes than spin-polarized calculations.Comparing the cell volumes of rutile VO 2 and V 2 O 5 , the additional oxygen and the concomitant transformation to the layered structure of the V 2 O 5 phase leads to an increase of the volume per formula unit by a factor of 1.5 to 1.7 Similar to the c-lattice vector, the calculated volume is either underestimated by 2 % for SCAN and SCAN+U, and by 4 % for SCAN+rVV or overestimated by 12 % for PBE and 13 % for PBE+U.
One of the fundamental properties of the electronic structure is the value of the band gap between the valence and the conduction band.The rutile VO 2 phase is known to be metallic, while V 2 O 5 is a semiconductor with an electronic band gap of about 2.3 eV [32].As shown in Table 1, all DFT functionals except the ferromagnetic solutions for PBE+U and SCAN+U correctly predict the VO 2 phase to be metallic.In the latter two cases, a 0.42 eV and a 1.04 eV band gap is opened for structurally relaxed cell geometries.Quite clearly the opening of a gap also depends on structural details since no gap is opened using PBE+U at U = 2 eV and experimental structural data [10].Interestingly, the gap in semiconducting V 2 O 5 is well described by all used functionals, with the largest deviation, 0.4 eV, found for SCAN+rVV.
In order to compare the performance of all functionals and spin configurations, two more key figures are displayed in Table 1: the magnetization energy ∆E FM-NM for rutile VO 2 and the reaction enthalphy for V 2 O 5 with respect both to nonmagnetic and ferromagnetic rutile VO 2 (H V2O5 f ).Table 1 shows that the calculated reaction enthalpy is always lower for the non spin-polarized calculations due to constraining all magnetic moments in the VO 2 phase to zero, which costs 79 meV to 720 meV per formula unit depending on the chosen DFT functional.The PBE functional yields the smallest value for ∆E FM-NM , and is increased by 136 meV for SCAN and by 242 meV for PBE+U.SCAN+U combines both effects and therefore leads to the largest value of 720 meV for ∆E FM-NM .As a result, non-magnetic calculations overestimate the enthalpy of reaction for V 2 O 5 by 520 meV to 790 meV.Employing spinpolarized PBE+U, SCAN and SCAN+rVV functionals reduces this error to 100 meV to 180 meV per V 2 O 5 unit.
The PBE functional overestimates H V2O5 f even for a ferromagnetic treatment, which is the consequence of two facts: First, the heat of formation is the lowest among the used functionals when performing non-magnetic calculations and secondly, the magnetization energy for the PBE functional is significantly lower compared to PBE+U and SCAN.On the other hand, the spin-polarized SCAN+U functional strongly reduces H V2O5 f as compared to SCAN and underestimates the experimental value by 890 meV.A similar decrease of H V2O5 f is also visible inspecting the values for the PBE and PBE+U functionals, where the PBE value is quite strongly decreased by taking into account on-site V-d coulomb interactions.
Investigating the electronic structure of metallic rutile VO 2 with the help of Wannier projections onto atomic-like vanadium d orbitals allows for a more explicit discussion.Figure 2 shows the total density of states (DOS) and projections to vanadium atomic orbitals calculated with the above set of functionals and spin configurations.Since the spin-polarized functionals, PBE, SCAN, SCAN+rVV yield similar projections as well as the PBE+U and SCAN+U we will only discuss the differences between the PBE and PBE+U functionals.Considering non-magnetic calculations, all functionals yield a band structure similar to the one obtained from PBE calculations (see lowest panel in the Figure 2).As shown, the band structure includes oxygen 2p bands in a range from −8 eV to −1.6 eV and vanadium 3d bands between −1.6 eV to 6 eV.Due to the crystal field the vanadium d bands are split into lower t 2g states ( d xz , d yz and d xy orbitals) and higher e g states ( d z 2 and d x 2 −y 2 orbitals).Including spin-polarization the splitting of the spin-up (↑) and spin-down (↓) channels, pushes the t 2g↓ states above the Fermi level.Consequently, all electrons in the vanadium 3d bands are located in t 2g↑ states.Including on-site coulomb correlations via an Hubbard U splits the t 2g↓ further and separates the d xy orbital from the d xz and d yz orbitals, and hence a 0.42 eV energy gap is opened.Here, all electrons in the V-3d band are found in the d xy orbital which coincides with the edge of the O-2p band.These results should be directly compared to state-of-the-art dynamical meanfield theory (DMFT) benchmark calculations as found e.g. in Fig. 1 of Ref. [13] (Note, d xy in the present work corresponds to a 1g in [13], d xz to e π g (1), d yz to e π g (2), and the remaining orbitals to e σ g ).A visual comparison directly reveals the improper description of the orbital occupancies for the spin-polarized PBE+U case, where a monoclinic-like  [13]).However, this spurious behaviour is only observed for larger values of U (U ≥ 2 eV) and using the relaxed optimized lattice where the c-direction comes out too large.
The SCAN+U functional describes the structural and electronic properties of the rutile VO 2 even worse: The rutile c-lattice vector is severely overestimated which leads to an even stronger distorted unit cell, the spurious band gap is increased to 1.04 eV, the respective orbital occupations are in disagreement with DMFT benchmark calculations and the calculated reaction enthalpy H V2O5 f deviates by more than 200 % from experimental value (1.3 eV).Therefore, the SCAN+U functional has been discarded in the subsequent surface studies.

VO 2 surface terminations
In the following our evaluation of the performance of the various functionals together in conjunction with the treatment of the spin configuration is extended to surfaces.The primary focus is on the (110) surface but to allow for a comparison to previous work which used (NM) PBE [15] or (FM) PBE+U [16] functionals, also these terminations have been considered briefly.The calculated surface energies for five different surface terminations for our chosen DFT functionals and spin configurations are shown in Table 2, including the already reported values for the PBE and PBE+U functionals.For the PBE functional, our values for the (100) and (001) surfaces are in perfect agreement with Mellan and coworkers [15], but the present work differs significantly, by up to 25 %, for the more open (110) surface.The origin of this difference can be traced to a differing computational setup, namely a different number of core electrons in the chosen PAW potential, and a different thickness of the slab.The surface free energies in [15] have been calculated using a four layer slabs only, whereas the slab thickness in our calculations was varied from five layers to eight layers to ensure a proper convergence of the surface energies.For the (NM) PBE+U functional our results are again in a perfect agreement with the values reported by Wahila [16], except for the (101) surface where our value is, surprisingly, larger by 90 %.This difference is presumably due to a different choice of the U value (U = 2 eV here and U = 3.25 eV in [16]), which according to our experience opens a band gap in the middle of a (101) slab.Our results imply that the calculated surface energies only depend slightly on the spin configuration, aside from two special cases.In the first case, values calculated with the PBE+U functional depend much more strongly on the spin configuration, which we blame on the spurious electronic instability in the V-t 2g bulk states for (FM) spin configurations, as discussed above.While the (NM) PBE+U functional yields an occupation of all the orbitals comprising the V-t 2g bands, the (FM) PBE+U functional almost exclusively occupies only the V-d xy orbital, as shown in Figure 2.However, using unrelaxed slabs close to the respective bulk geometries surface energies are again much less affected by the spin configuration, resulting in only an 13 % decrease for the (NM) configuration, which confirms that the surface free energy differences are indeed caused by the spurious electronic transition (gap opening) in the bulk and respectively optimized geometries.In the second case, (NM) non spin-polarized calculations show a systematic ∼20 % decrease of the (001) surface free energy for all the functionals.This decrease is not observed for unrelaxed slabs, where the calculated surface energies are more or less the same for (NM) and (FM) spin configurations.The additional stabilization of the relaxed slabs is accompanied by a change in the surface V-t 2g states depicted in Figures 3(a-b).As shown in panel 3(a), (NM) non spin-polarized calculations show a peak centered around −0.7 eV below E F formed almost exclusively by V-3d states located close to or at the surface.This is not the case for (FM) spin-polarized calculations where the V-3d DOS is quite similar for both bulk and at the surface vanadium atoms as easily detectable in 3(b)).As the (001) surface is the only surface termination which comprises parallel chains of alternating VO 2 octahedra and tetrahedra (see panel 3(c)), we suggest that the change in the electronic structure is driven by this specific surface structure.
In summary, the values of the surface free energies shown in Table 2 confirm that a VO 2 (110) surface termination is the most stable one, independent of the chosen DFT functional and spin configuration.PBE and (NM) PBE+U functionals yield rather low values, which is in line with a previous conclusion that PBE underestimates the surface energy for correlated materials [41].But all other tested functionals provide more reasonable values.Concluding this paragraph, we note that a comparison of surface free energies unfortunately reveals a strong dependence on the chosen functional, which must be kept in mind especially when discussing absolute numbers.We also suppose that similar effects might be present in other surface orientations not treated in the present work.bare (110) surface, resulting in the (2 × 1) buckled superstructure shown in Figure 4.The surface layer shows an unpaired zig-zag line pattern, similar to the monoclinic M2 phase [7].The relative height difference in this buckled row depends on the respective functional and ranges from 0.30 Å to 0.43 Å.An instability towards buckling can already be observed for the bulk phases of VO 2 : On the one hand, a V-V pairing is present along the [001] direction of the monoclinic ground state.One the other hand, several DFT functionals, such as PBE, predict an imaginary DFT phonon mode for the rutile structure [42], where the oxygen atoms are shifted along the [001] direction with respect to the vanadium atoms.To exclude the possibility that the buckling at the surface is driven by this artificial bulk mode, we recalculated the structure where the bottom layer was fixed at its bulk structure and the remaining atoms were not allowed to relax in the [001] and [1][2][3][4][5][6][7][8][9][10] direction.Even under these constraints, the buckled surface has again a lower energy: buckling stabilizes the surface by 40 meV, 15 meV, 9 meV and 11 meV per (2 × 1) slab using the (NM) PBE, PBE+U, SCAN and SCAN+rVV functionals, respectively.These results reinforce the above assumption that this buckling is indeed a surface effect and not driven by instabilities in the bulk phase.Fully relaxing the slabs increases the stability induced by the buckling by 72 meV, 39 meV, 45 meV and 59 meV per the (2 × 1) slab, respectively, for the functionals above.The buckling is accompanied by a change in the projected DOS onto V-3d orbitals of the upper vanadium atoms in the buckled row.The occupation of the projected V-3d DOS is decreased by ∼0.3 electrons (PBE) as compared to an unbuckled surface.To the opposite, spin-polarized calculations with (FM) spin ordering reveal a different trend: while PBE, SCAN and SCAN+rVV functionals strongly suppress the stability of the buckled surface  and transform it back to a bare (110) surface, the buckling is recovered at the PBE+U level, with a 26 meV energy gain for fully relaxed slabs.Since the small energy gain due to the buckling might render an experimental verification difficult, we performed a molecular dynamics (MD) simulation at experimental temperatures to evaluate the thermal stability of the buckling.The simulation for 100 ps was performed employing the non-magnetic PBE functional, which shows the largest energy difference between the buckled and the bare (110) surface, a temperature of 350 K and a time step of 1 fs.We monitored the most prominent feature of the buckled surface, namely possible switches in the z-positions (i.e.along the surface normal) of the V-V pairs.In order to decrease the computational cost we reduced the thickness of the slabs to four layers where the bottom layer was fixed and the k-points grid was reduced to 2 × 2 × 1 points.The results of the MD simulations are shown in Figure 4 as the number of switches at the y axis and the buckling size at the x axis.The plot also shows a fit of the MD data (black line) with three Gaussian functions, related to the buckled (red line), non-buckled (green line) and flipped buckled (blue line) configurations.The sum of the three Gaussian fits (gray line) shows that the buckled configuration is thermally smeared since there are no maxima at the positions of the Gaussian functions related to the buckled configurations.Furthermore, the average flipping time from the up-down to downup configurations was calculated to be 260 fs.Thus, we conclude that the time dependent stability of the buckling is much shorter than the resolution limit of atomically resolved STM experiments, as imposed by the typical scanning rate.
To shed some more light on the origins for the buckling we investigated the influence of the buckling on the electronic structure with the help of Wannier projections.Figure 5 shows the results using the SCAN functional.The left, middle and right panels show projections onto surface vanadium atoms at unbuckled, buckled-up and buckled-down sites.An inspection of these projections reveals that a V-d yz orbital in the buckled surface is less occupied as compared to to the unbuckled surface and the occupation of a V-d xz orbital depends on the z-position (up or dn) of the buckled vanadium atoms.The Wannier projections reveal that the impact of the surface buckling on the electronic structure is qualitatively different from  changes induced by MIT.As it was pointed out by Eyert [7], formation of vanadium pairs in the monoclinic phase is accompanied by an increase in the occupation number of the d xy orbital.Note that in the present work the notation of d xy and d x 2 −y 2 is interchanged.However, the buckled and the unbuckled surface show an opposite behavior: As shown in Figure 5, the occupation number of the d xy orbital localized at both buckled and unbuckled sites is significantly reduced compared to the rutile bulk (see Figure 2, top panel).
Since the occupation of the V-3d states is quite sensitive to buckling, we also looked into to what extent a varying oxygen coverage influences the size of the buckling, comparing bare (clean), 25 %, 50 % and 100 % covered VO 2 surfaces as depicted in Figures 6(a-d).Figure 7 clearly shows that all non-magnetic functionals predict a stable buckled surface and that its stability increases with oxygen coverage.The energy gain due to the buckling ranges from 39 meV (PBE+U) to 72 meV (PBE) for the bare surface, as already discussed above, and increase with oxygen coverage, adding further stability which ranges from 61 meV (PBE+U, SCAN) to 83 meV (SCAN+rVV) per the (2 × 1) slab for a fully-covered surface.Quite evidently van der Waals interactions included in the SCAN+rVV functional have only a negligible impact on the stability of the surface buckling.Since all spin-polarized functionals except PBE+U predict an unbuckled surface, these negative values are omitted in Figure 7, and only the (FM) PBE+U value of 26 meV is included.However, the trend changes with increasing the oxygen coverage to 50 %, which compared to non-magnetic runs yields a further 20 % (PBE), 60 % (PBE+U), 70 % (SCAN), and 40 % (SCAN+rVV) stabilization of the surface buckling, which ranges from 100 meV to 140 meV per the (2 × 1) slab on an absolute scale.Increasing the oxygen coverage to the fully-covered limit, the buckled surface is most preferred for the (FM) PBE+U functional which adds a 510 meV extra stability per the (2 × 1) slab.The other functionals, PBE, SCAN and SCAN+rVV also prefer a buckled surface by 104 meV to 137 meV per the (2 × 1) slab without any significant difference between the non-magnetic and ferromagnetic spin configurations.We attribute the comparably large energy gain for the (FM) PBE+U functional to the opening of the gap in the bulk and hence to a reduced screening of the induced charge imbalances upon oxygen adsorption.Hence, this buckling might play a non-negligible role for the stabilization of oxygen-rich surface reconstructions proposed in recent experimental studies [19,16] and even more on surfaces of the insulating M phase of VO 2 .
Since the experimentally observed (2 × 2) superstructure [18] cannot be explained by the surface buckling only, we investigated the stability of a large set of surface reconstructions obtained by optimizing preselected random structures for a number of surface stoichiometries and each employed functional and spin configuration.The detailed procedure is outlined in the Supplementary Information.In the following, a given surface stoichiometry is defined by the number of vanadium and oxygen atoms added on top of bare rutile VO 2 (110) surface.Independent of the chosen functional and spin configuration, the five most stable terminations which resulted from optimizing the preselected random structures are depicted in Figures 8(bf).All of them are based on a general pattern shown in   [43].This reconstructed surface superstructure consists of alternating octagonal and hexagonal rings made of VO 2+x polyhedra as most easily visible in the top view of panel 8(a).However, the surface free energy that results from this termination is too high so that a further oxidation is required to stabilize this superstructure by adding oxygen atoms in between the octahedral rows of the subsurface layer (see the blue circle).The first stable (2 × 2) ring superstructure derived from the general pattern is shown in Figure 8(b).The modification proceeds in two steps: First, the bridging vanadium tetrahedra (see the green circles) are removed and second, two additional oxygen atoms are adsorbed on the ring pattern, which pushes the vanadium tetrahedra towards the vacuum as indicated by the the blue arrow.These modifications lead to the formation of a stable, disconnected ring structure with a V 4 O 13 surface stoichiometry and showing only VO x tetrahedra.This ring termination is quite flexible and can be easily shifted by a quarter of the c-bulk lattice vector as indicated by the black arrow in Figure 8(c).As a result, the lower VO x tetrahedra are additionally bound to a second oxygen ad-atom which changes their coordination geometries to square pyramids.This change is most easily visible by comparing the features within the black circles in Figures 8(b) and  (c).Both disconnected ring structures in Figures 8(b,c) lead to a rectangular STM pattern (see Figure S2(a) in the Supplementary information) where the bright features are due to the topmost oxygen atoms.The rel-ative stability of these ring structures depends on the choice of the functional and will be discussed further below.Proceeding with the ring structures, another variant is shown in Figure 8(d), which descends from the general pattern in the following way: First, the ring structure is shifted by a quarter of the [1 10] bulk lattice vector (see the red arrow in the top view) so that the bridging tetrahedra are aligned with upright octahedra of the subsurface layer as marked by the orange circle.Second, this shift of the ring termination relative to the subsurface layer is accompanied by the oxidation of different vanadium tetrahedra.All together this leads to a termination with an overall V 5 O 15 surface stoichiometry displaying a zig-zag STM pattern (see Figure S2(b) in the Supplementary information).The next stable superstructure shown in Figure 8(e) is generated by shifting the previous structure by an additional quarter of the [001] bulk lattice vector (see orange arrow), which again causes the bridging tetrahedra to bind to two subsurface octahedra, as easily seen by comparing the features inside the orange circles in Figures 8(d,e).However, here, one oxygen atom is removed changing the overall surface stoichiometry to V 5 O 14 and causing the tetrahedral coordination geometry not to switch to a square pyramid as has been the case for the V 4 O 13 ring structures.
Since the ring terminations are supported on a bare surface, also oxygen atoms might be missing in the connecting subsurface layer which changes the local coordination of the affected vanadium atoms from octahedra to square pyramids.This is shown in Figure 8(f) where we consider a ring structure similar to the shifted V 4 O 13 ring (see Figure 8(c)) where the hexagonal rings are connected by bridging tetrahedra.The side view shows the stepwise removal of subsurface oxygen atoms surrounding two vanadium sites (see the green and violet circles in Figure 8(f)), which leads to a V 5 O 13 and a V 5 O 12 surface stoichiometry, respectively.
Our results reveal that the relative stability of the two V 4 O 13 surface terminations shown in Figures 8(b,c) depends on the chosen functional.The energetic differences shown in Table 3 clearly indicate that the (FM) spin-polarized PBE and PBE+U functionals prefer the purely tetrahedral ring structure (b) by 114 and 385 meV per (2 × 2) supercell while all the other functionals find the shifted modification (c) as more stable by 4 meV to 324 meV.
In order to compare the stability of the superstructures with different surface stoichiometries, we plot the surface free energy of the most stable terminations in Figures 9 and 10 for ferromagnetic (FM) and nonmagnetic (NM) calculations, respectively.Furthermore, we compare the ring terminations with a bare (110) surface (black line) and simple oxygen ad-atom (110) surface models shown in Figure 6.The green '1O',  Oxygen adsorption on a rutile (110) surface has been already studied by Mellan et al. [15] performing non spin-polarized PBE calculations.In their work, the Γ = +1 and Γ = +2 surfaces are identical to our 50 % and fully covered VO 2 (110) terminations depicted in Figures 6(c,d).However, our calculated adsorption energies for both coverages are ∼0.46 eV larger than the corresponding (i.e. the uncorrected) values for the Γ = +1 and Γ = +2 surfaces (see Figure 8 in ref. [15]).This rather large discrepancy in adsorption energy is caused again mainly by a different computational setup.The use of differing PAW potentials causes a 0.26 eV difference and the different slab thickness increasing this value by another 0.07 eV.Furthermore, the models in the present work are additionally stabilized by the surface buckling, which has been neglected in the work of Mellan et al. [15] which increases the difference by ∼0.06 eV for both coverages.The remaining small discrepancy is attributed to other effects emerging from the different computational setup.
Figures 9 and 10 show the calculated surface free energies of the bare VO 2 (110) surface (black lines), the oxygen ad-atom phases with different oxygen coverage (green lines) and the ring terminations depicted in Figures 8(b-f) for varying oxygen chemical potential.Concerning the V 4 O 13 ring structure only the surface free energy of the most stable variant is included.Independent of the chosen spin configuration all graphs show that all the functionals considered determine the ring terminations to be more stable than (PBE, PBE+U, SCAN) or at least similarly stable (SCAN+rVV) as the oxygen ad-atom phases over the wide range of chemical potentials.However, the relative stability of the ring terminations with respect to the oxygen ad-atom phases is functional dependent.To quantify this dependency we consider the difference between the crossing points of both the fully covered oxygen ad-atom (110) surface (light green line) and the V 5 O 14 ring phase (cyan line) with the bare stoichiometric (110) surface, which determine the effective adsorption energy of the additional oxygen atoms that are incorporated in the superstructure.Independent of the spin treatment, the PBE and PBE+U functionals determines this difference to fall between 0.4 eV to 0.5 eV, SCAN prefers the V 5 O 14 ring by 0.1 eV and the SCAN+rVV functional renders the difference in crossing points to below 0.01 eV.The graphs also show that the (FM) PBE, (NM) PBE and (NM) PBE+U functionals prefer the ring structures over the whole range of the oxygen chemical potentials while the other functionals and spin configurations yield an energy window in the low oxygen chemical potential region where the oxygen ad-atom phases would be preferred.On the other hand, the (FM) SCAN, (FM) SCAN+rVV and (NM) SCAN+rVV functionals show a preference for the ring terminations only in regions where the VO 2 bulk phase is calculated to be thermodynamically unstable.
The strong stabilization of the ring structures for the PBE and PBE+U functionals also hints towards an existence of additional polyhedral surface terminations that are more stable than the ad-atom phases.Regarding the V 4 O 13 and V 5 O 15 rings, most of the functionals do not exhibit a significantly different stability values except for the (NM) SCAN, SCAN+rVV and (FM) PBE.While the V 5 O 15 ring is more stable by 50 meV (SCAN) and 80 meV (SCAN+rVV) than the disconnected ring, the (FM) PBE functional prefers the connected V 5 O 15 ring by 50 meV.
For the ring structure on the reduced oxygen missing subsurface layer (panel 8(f) we find again a strong dependency on the chosen functional.The graphs 9 and 10 find no additional stability at the SCAN or SCAN+rVV level, while these structures are preferred at the PBE and PBE+U level under strongly reducing conditions, see yellow dashed lines.Particularly the (NM) PBE+U functional prefers this structure for oxygen chemical potentials in the range of −1.5 eV to −2.8 eV.A enhanced stability is also found at the (NM) and (FM) PBE level, showing that a doubly reduced subsurface layer would be stable for an oxygen chemical potential below −1.9 eV and −1.8 eV respectively.A stability window ranging from −1.23 eV to −1.89 eV is also found for (FM) PBE+U setups.The qualitative agreement between the PBE and PBE+U functionals can be understood by considering the existence of an electronic gap for states connected with the surface layer which will be discussed in the following section.The other functionals determine these reduced subsurface layer structures to be unstable with respect to both the other other ring structures and ad-atom phases.Comparing the reduced V 5 O 13 and V 5 O 12 ring terminations reveal yet another functional dependent property.Considering their crossing points with the bare (110) surface line, the PBE and PBE+U functionals show crossing to occur at a lower oxygen chemical potential for the more reduced ring termination with a V 5 O 12 surface stoichiometry.However, these crossing points are found either roughly at the same oxygen chemical potential (SCAN) or they are even reversed (SCAN+rVV), revealing a counter-intuitive opposite trend with respect to the ad-atom phases, see green lines in Figures 9 and 10.Vacancies at PBE and PBE+U level stabilize the surface which is not the case for the SCAN functional.
One of the most striking differences between the spin-polarized and non spin-polarized calculations for the surface free energies concerns the absolute values of the effective adsorption energies.This difference can be illustrated best for the ring termination with a V 5 O 14 surface stoichiometry (Figure 8(e).While the non-magnetic calculations yield effective adsorption energies, depending on the chosen functional, between −1.88 eV to −2.06 eV, (FM) spin-polarized calculations shift the resulting values by 0.20 eV, 0.60 eV, 0.61 eV and 0.57 eV to the right (less negative values) using PBE, PBE+U, SCAN or SCAN+rVV.Note that similar shifts (±0.1 eV) are also observed for the ad-atom phases (green lines in Figures 9 and 10).This energy shift can be attributed to higher energies for bulk rutile VO 2 when performing non-spin polarized calculations which are due to the neglect of local magnetic moments, see Table 1.
Considering a possible opening of a gap in the surface DOS of the oxygen rich ring terminations it is tempting to inspect the relation between the V 2 O 5 bulk phase and the ring terminations more closely.In Figure 11 we compare the projected density of states of bulk VO 2 and V 2 O 5 phases with the oxygen ad-atom phase at full coverage (Figure 6(d)) and the V 4 O 13 ring.major difference between the electronic structure of VO 2 and V 2 O 5 bulk phases is found for the occupation of the V-3d states.While the Fermi level cuts the V-3d bands in the VO 2 phase, the V 2 O 5 phase shows an electronic band gap of about 2.0 eV between the O-2p and V-3d states.As can been seen in Figure 11, the projected surface DOS of the ring structure resembles rather closely the DOS of the V 2 O 5 bulk phase, revealing an unoccupied V-3d band and hence this surface termination cannot form local magnetic moments.Figure 11 also shows that the terminating ring layer is insulating with a 1.9 eV band gap, which is in a good agreement with the 2.3 eV obtained for the V 2 O 5 bulk (see Table 1).It should also be noted that the surface free energies of the ring terminations correlate with the calculated oxidation enthalpy of bulk VO 2 .As visible from the surface free energy graphs above, the difference between the effective adsorption energy of the V 5 O 15 ring phase and the reaction enthalpy for a V 2 O 5 formation by oxidizing VO 2 is below 0.1 eV in all cases, while the absolute errors in calculated oxidation enthalpies range from 0.1 eV to 0.79 eV.These facts indicate, that the ring terminations are related to a monolayer of V 2 O 5 (001) in agreement with our previous work [18].This finding also implies that a correct evaluation of the stability for both phases is important to properly describe the energetics of the oxygen rich surface terminations, a criterion which seems to be matched best by spin-polarized PBE+U, SCAN and SCAN+rVV calculations.However, the application of the (FM) PBE+U functional yields wrong orbital occupations and strongly overestimates the c/a ratio of the bulk rutile VO 2 phase.From this perspective the spin-polarized SCAN and SCAN+rVV functionals seem to be a better choice for performing calculations on VO 2 surfaces.

Summary
In this work we studied the performance of spin polarized and non-spin polarized GGA, GGA+U meta-GGA and meta-GGA+U functionals used to characterize bulk properties, stoichiometric surface terminations and oxygen-rich (2 × 2) (110) reconstructions of the rutile VO 2 phase.We find that surface energies of particular surface orientations depend on the spin treatment due to complex electronic structure of both the surface states or the bulk states located in t 2g orbitals.Considering a stochiometric (110) surface, we identified a ground state with a buckled surface layer, resulting in a (2 × 1) superstructure.The relative stability of the buckling depends on the chosen functional and becomes larger upon increasing the oxygen coverage which in turn reduces the occupation of the V3d band near the surface.At high oxygen coverage, the (110) surfaces are terminated with V 2 O 5 related ring structures.We identified three important factors that influence the calculated stability of the these reconstructions.First, the PBE and PBE+U functionals prefer quite open terminations loosely based on tetrahedral VO 2 polyhedra.This preference is already visible in the bulk calculations of the vanadium pentoxide phase, where the layer spacing is overestimated while the SCAN and SCAN+rVV functionals underestimate it (see Table 1).The leads to a large number of polyhedral superstructures that are stable at the PBE level, but unstable with the more advance SCAN meta-GGA functional.In addition, our results show that the V 5 O 15 ring termination is closely related to the vanadium pentoxide, which indicates that a correct description of the relative stabilities of the rutile VO 2 and V 2 O 5 phases plays an important role for the calculation of the surface free energies.Our data show that spin-polarized calculations are superior to the non-spin polarized ones for the oxygen covered surfaces, yielding in all cases (except the SCAN+U) the stability of the V 2 O 5 phase with respect to VO 2 closer to the experimental findings.The third factor that modifies the surface free energy is related to the surface buckling.Yet, while this buckling can change the adsorption energy by as much as 0.5 eV for a (2×1) surface cell when using PBE+U Hubbard corrections for a (FM) spin configuration, this effect is much less pronounced for the other functionals.In summary, a ring structure reconstruction of the (110) surface is either energetically preferred or at least as stable as the adsorption phase at high oxygen chemical potential, independent of the chosen functional and spin treatment.
First principles studies of the electronic and structural properties of the rutile VO 2 (110) surface and its oxygen-rich terminations After that, all of them were pre-optimized, and the 10 % most stable ones fully relaxed, which led to a set of 100 structures, similar to the one from panel (b).At the next step, an additional substrate layer was added as shown in panel (c).Since two substrate layers were too thin according to our convergence studies, the last step consisted of adding two more substrate layers, as shown in panel (d).Concluding, a full and accurate relaxation with all functionals using both spin-polarized (FM) and non-spin polarized (NM) electronic configurations was performed.
we considered V 8 O 18 and V 8 O 20 surface stoichiometries, which ist the the same as for the whole VO 2 trilayer with 2 and 4 additional oxygen atoms, respectively.With our approach we obtained the fully and half oxygen covered (110) surface respectively and therefore concluded that the configuration space is sampled sufficiently well, even though there is no guarantee for finding at the given stoichiometry the most stable structure.
In our work we first considered strongly oxidized surface stoichiometries V n O 2n+2 to V n O 2n+5 for n ranged from 4 to 8 which corresponds to the more than a half of an additional VO 2 trilayer.Nevertheless, in the case of surface terminations with n ≥ 6, we either observed formation of the second surface layer or the formation of a (110) termination with a V vacancy.Both these superstructures turned out to be unstable with respect to the bare (110) surface and other ring-like terminations.Properties of the rutile VO 2 (110) surface and its oxygen-rich terminations

Figure 1 .
Figure 1.The structures of V 2 O 5 (left panel) and rutile VO 2 (right panel).The dashed lines in the left panel mark the inter-layer van der Waals bonds along the c-direction between V 2 O 5 layers.The yellow arrows in the right panel denote the z -directions of the individual VO 6 octahedra's local coordinate systems, which are mutually orthogonal to the c-direction of the rutile VO 2 unit cell.

Figure 2 .
Figure 2. Density of states (DOS) of the rutile (R) VO 2 phase calculated with different DFT functionals and projected onto Wannier orbitals localized on vanadium atoms.

Figure 3 .
Figure 3. Layer-resolved V-t 2g DOS of a 7-layer VO 2 (001) slab and employing PBE.Surface denotes the DOS projected onto V atoms located both in the surface and subsurface layers, while bulk denotes the contributions of the three inner bulk-like layers.Panels (a) and (b) show the respective DOS for NM and FM spin configurations.The inset (c) shows the (001) surface structure comprising parallel chains of alternating VO 2 octahedra and tetrahedra.

Figure 4 .
Figure 4. Left side: Buckled (2 × 1) reconstruction obtained from simulated annealing for non spin-polarized calculations.The size of the buckling in the surface layer amounts to 0.2 Å, and is accompanied by a decrease of the occupation (0.3 electrons) of the V3d states of the upper vanadium atoms indicated by the black arrow.Right side: Distribution of height differences between the V-V pairs during a molecular dynamics simulation consisting of 10 5 steps with a time step of 1 fs performed at 80 °C (black line) including a fit with three Gaussian functions (red, blue and green lines) corresponding to the up-down, non-buckled and down-up configurations of the V-V pairs, respectively.The gray line is the sum of all fit functions.

Figure 5 .Figure 6 .
Figure 5. (NM) non spin-polarized SCAN DOS projected onto Wannier functions localized at surface vanadium atoms of the bare and buckled (110) surface.Only the V-3d contributions of the surface V atoms of the unbuckled (Bare(110)) and upper V (Buckled (110) up) respectively lower V (Buckled (110) dn) of the buckled VO 2 surface are shown.

Figure 8 (
Figure8(a) which displays a V 5 O 12 surface stoichiometry and is structurally similar to the (3 × 1) surface reconstruction of SrTiO 3 (110)[43].This reconstructed surface superstructure consists of alternating octagonal and hexagonal rings made of VO 2+x polyhedra as most easily visible in the top view of panel 8(a).However, the surface free energy that results from this termination is too high so that a further oxidation is required to stabilize this superstructure by adding oxygen atoms in between the octahedral rows of the subsurface layer (see the blue circle).The first stable (2 × 2) ring superstructure derived from the general pattern is shown in Figure8(b).The modification proceeds in two steps: First, the bridging vanadium tetrahedra (see the green circles) are removed and second, two additional oxygen atoms are adsorbed on the ring pattern, which pushes the vanadium tetrahedra towards the vacuum as indicated by the the blue arrow.These modifications lead to the formation of a stable, disconnected ring structure with a V 4 O 13 surface stoichiometry and showing only VO x tetrahedra.This ring termination is quite flexible and can be easily shifted by a quarter of the c-bulk lattice vector as indicated by the black arrow in Figure8(c).As a result, the lower VO x tetrahedra are additionally bound to a second oxygen ad-atom which changes their coordination geometries to square pyramids.This change is most easily visible by comparing the features within the black circles in Figures8(b) and (c).Both disconnected ring structures in Figures8(b,c) lead to a rectangular STM pattern (see FigureS2(a) in the Supplementary information) where the bright features are due to the topmost oxygen atoms.The rel-

Figure 8 .
Figure 8. Perspective and top views of the resulting (2 × 2) superstructures obtained from an optimization of random structures with different surface stoichiometries.The structures descend from an underlying pattern shown in panel (a) which is structurally similar to the (3 × 1) reconstruction of the SrTiO 3 (110) surface [43].

Figure 9 .V 4 O 13 PBETable 3 .
Figure 9. Calculated surface free energies as a function of the oxygen chemical potential for various (2 × 2) VO 2 (110) surface terminations and a FM spin configuration.

Figure 10 .
Figure 10.Calculated surface free energies as a function of the oxygen chemical potential for various (2 × 2) VO 2 (110) surface terminations and a NM spin configuration.

Figure 11 .
Figure 11.Vanadium and oxygen projected DOS of the VO 2 , V 2 O 5 phases, and projections to surface atoms of a VO 2 (110) surface fully covered with oxygen ad-atoms and of the V 4 O 13 ring termination.The projected DOS was calculated with the non-magnetic PBE functional and is both normalized to a single atom and aligned at the upper edge of the O-2p band.

Figure S1 .
Figure S1.Optimizaton process of the random structures.At the beginning, a set of 1000 random structures was generated with the USPEX[1]  package as it is shown in panel (a).After that, all of them were pre-optimized, and the 10 % most stable ones fully relaxed, which led to a set of 100 structures, similar to the one from panel (b).At the next step, an additional substrate layer was added as shown in panel (c).Since two substrate layers were too thin according to our convergence studies, the last step consisted of adding two more substrate layers, as shown in panel (d).Concluding, a full and accurate relaxation with all functionals using both spin-polarized (FM) and non-spin polarized (NM) electronic configurations was performed.

5 1. 3 .z 2 d xz d yz d xy d x 2 −y 2 PBEFigure S3 .
Figure S3.Projected density of states onto Wannier orbitals localized at vanadium atoms in the rutile (R) phase, calculated with several DFT functionals using relaxed (two left columns) and experimental structures (two right columns).0 marks the Fermi energy, t 2g states are formed by d xz , d yz and d xy orbitals, d z 2 and d x 2 −y 2 form the e g states.

Table 1 .
Structural and energetic details of the rutile VO 2 and V 2 O 5phases, calculated with several DFT functionals for ferromagnetic (FM) and non-magnetic (NM) spin configurations at respectively optimized relaxed structures.c/astandsfor the ratio of the c and a lattice parameters and V for the volume of the unit cell (in Å3 ) .Eg denotes the band gap (in eV) and ∆E FM-NM the magnetization energy (in meV) defined as energy difference between the FM and NM spin configuration.H V 2 O 5f is the reaction enthalpy for forming the V 2 O 5 phase with respect to rutile VO 2 .The volumes, ∆E FM-NM and H V 2 O 5 f are given per formula unit.by 0.42 eV from the d xz and d yz orbitals, rendering rutile VO 2 a semiconductor.This clearly contradicts the benchmark DMFT results for rutile VO 2 (R) where the DOS does not show this splitting and the system remains metallic (see upper panel (a) of Fig. 1 in

Table 2 .
Calculated surface free energies in meV/ Å2 for various VO 2 surface orientations.