Determining Hubbard U of VO2 by the quasi-harmonic approximation

Longjuan Kong(孔龙娟), Yuhang Lu(陆雨航), Xinying Zhuang(庄新莹),Zhiyong Zhou(周志勇), and Zhenpeng Hu(胡振芃)

School of Physics,Nankai University,Tianjin 300071,China

Keywords: quasi-harmonic approximation,vanadium dioxide,first-principles calculation,Hubbard U

1.Introduction

Strongly correlated transition metal oxides exhibit a wide range of intriguing physical properties, including magnetoresistance,superconductivity,and metal–insulator transition(MIT),which arise from their strong electron–electron correlation effects.[1–6]These materials have attracted significant attention in recent years and hold great potential for various applications.[7–11]Among them, vanadium dioxide VO2stands out as a representative material that has been extensively studied due to its remarkable feature known as the metal-to-insulator transition.[12]This transition was first discovered by Morin in 1959 and occurs around 340 K, accompanied by a simultaneous structural change from the high-temperature conducting rutile (R) phase to the lowtemperature insulating monoclinic(M)phase,with a band gap of approximately 0.6 eV.[12–14]The MIT of VO2near room temperature is associated with rapid changes in its structural,electrical,optical,thermal,and magnetic properties.[15–19]The reversible nature of this transition, along with the relative ease of synthesis,[20–23]makes VO2an ideal material for various applications, including smart windows,[24,25]photoelectric switches,[26]information storage,[27]and thermal control of spacecraft.[28]Moreover,the unique MIT behavior of VO2has gained renewed attention with the development of new ultrafast microscopy techniques.[29,30]The MIT in VO2,with its drastically changed physical properties, can be triggered by different stimuli such as temperature, pressure, light, electric field,and doping.[31,32]

Furthermore,VO2not only provides a playground for exploring the peculiar MIT and corresponding applications, but paves the way for tackling strongly correlated systems.Unfortunately, the precise elucidation of the physical properties of materials with strongly correlated electronic behavior has continued to be a challenging subject in condensed matter physics.An inexpensive and reliable way to evaluate the physical features of with strong correlated d/f electrons is the calculations based on or beyond the density functional theory(DFT).As a consequence,extensive theoretical works have been performed, including but not limited to classical DFT within the local density approximation(LDA),[33,34]and the generalized gradient approximation (GGA),[35]hybrid functions,[36]GW(G refers to Green’s function and W denotes the dynamically screened Coulomb interaction),[37]and dynamical mean field theory (DMFT).[38]It is noticed that LDA and GGA cannot describe the electronic structure of M-VO2correctly.The GW approximation and the random phase approximation (RPA)have been proved rather useful to describe the excitation of materials.The DMFT may correctly describe the electronic properties of both R and M phases of VO2.Unfortunately,these three methods fall short for the calculations of forces and suffer from the high computational cost.Though hybrid functional may be used to calculate forces,it still suffers from the high computational cost and fails in describing the metallic behavior of R-VO2.[39]Generally,the DFT+Umethod is tuned to be extensively used to study ground state properties such as structural, electronic and thermal properties with an economic computational cost.The strong electron–electron Coulomb interaction in the 3d orbital of VO2and the drastic change in electrical conductivity following MIT make this system can be studied using Hubbard model with the consideration ofUparameter.[40]Taking HubbardUas an adjustable parameter, DFT+Umethod is easy to reproduce the nature of VO2,and it is suitable to calculate the total energy,atomic force,and stress accurately.These physical quantities,including the phase transition,mechanical properties,characteristics of surface and interface,and chemical reactions,are essential for better understanding of the richness of VO2.Usually, the conventional DFT, adding a HubbardUparameter representing the strong electron–electron interaction of strongly correlated materials, can cure the deficiencies of local or semilocal exchange-correlation functions and maintain relatively low computational cost.[41–43]

As was reported, the accuracy of outcomes by DFT+Uis tend to depend on the magnitude ofUparameter.In the need of the correct description of strongly correlated system properties, it is necessary to determine a reliableUvalue firstly.Hence, many methods have been proposed to determine the optimalUvalue.The widely used approach reported in literature for choosing effectiveUis adjusting calculated results to fitting experimental data(for instance,fitting the experimental bandgaps, seeback effect or electrical conductivity).[42,44,45]On the other hand, several pure theoretical perspectives have also been developed to determine theUnon-empirically, for example, the linear response method based on constrained DFT (CDFT) by Cococcioni and de Gironcoli1,[46]the constrained random-phase approximation(cRPA)method,[47]and self-consistent methods from the local electron density based on the Thomas–Fermi screening mode.[48]However, these methods have shortcomings more or less.For example, the limitation of determining HubbardUby fitting the experimental band gap data may lead to a hugeUin some materials,which is a lack of physical meaning.Meanwhile,several pure theoretical methods may either suffer significant expensive computation and fall short for Hellmann–Feynman force calculations,[37,49]or fail to obtain the effectiveUvalue correctly describing two VO2phases simultaneously.[50]Considering the efficiency and accuracy, the effectiveUvalue can also be indirectly appreciated by the comparing calculated phase transition temperatureTcbased on quasi-harmonic approximation (QHA) with that of experiment, from the thermodynamical aspect during MIT process.Actually already in 2014,Budaiet al.provided a full thermodynamic understanding the nature of VO2MIT by using x-ray and neutron scattering margins,revealing the importance of anharmonic phonons of lattice vibrations in phase-transition process.[51]The QHA,which takes into part of the anharmonic effects,may be a much more straightforward and computationally cheap method used to investigate thermodynamics properties, including the thermal expansion coefficients, the heat capacities, vibration entropy,and free energy.

In this paper,we propose a method to determine the optimal HubbardUvalue of VO2based on the transition temperatureTcusing the QHA from a thermodynamic perspective.We demonstrate that the obtainedUvalue is highly effective in capturing the correct physical properties of both the monoclinic (M) and rutile (R) structures of VO2.Focusing on the metal-to-insulator transition of vanadium dioxide, we choose to fit the phase transition temperature to determine the HubbardUin DFT+Ucalculations with the quasi-harmonic approximation.DifferentUvalues have been evaluated for the properties of either monoclinic or rutile phase,including structural information, band structure, as well as the phase transition temperature.After considering the electronic part of free energy, theUof 1.5 eV is proved to be an optimal one for Perdew–Burke–Ernzerhof(PBE)+Ucalculations of VO2.It can well reproduce the essential properties of VO2, which is efficient in the calculations.Overall, our method provides a reliable approach for determining the appropriateUvalue in DFT+Ucalculations of strongly correlated materials like VO2.

2.Theoretical methods and computational details

In this work, we performed the geometrical optimization and electronic properties’ calculations under the framework of density functional theory, as implemented in the Viennaab initiosimulation package (VASP).[52,53]The projector augmented-wave (PAW)[54]pseduo potential was employed to describe the interaction between the core and valence electrons.The exchange–correlation interaction was described by the PBE functional.[55]The electron wavefunction was expanded in a plane-wave basis set with an energy cutoff of 520 eV.In order to ensure a reliable energy comparison between two VO2phases, a reciprocal-grid density ofwas used, wherea,b, andcare the lattice constants in units of angstrom.Then, thek-point grids of 9×11×9 and 12×12×20 were determined for M and R cells,respectively.Since the previous studies have shown that the physical properties of strongly correlated system depend significantly on the magnitude ofU, we rationally tested different HubbardUvalues.The convergence criteria for force and energy are 10−3eV/˚A and 10−8eV,respectively.

The PHONOPY software[56,57]was employed to calculate the phonon related information at harmonic and quasiharmonic levels, including phonon dispersion, vibration free energy, the heat capacity, Helmholtz free energy, and the Gibbs free energy.The force constants of the supercell including 48 atoms for either M or R phase were calculated through the density-functional perturbation theory(DFPT).[58]A mesh grid of 40×40×40q-points was generated to make calculation more accurate.

The finite temperature-dependent physical quantities were obtained by using the QHA,[56,59]including vibration free energy,the heat capacity,Helmholtz free energy,and the Gibbs free energy.The temperature was assumed to indirectly affect phonon frequencies through thermal expansion and temperature dependence of phonon frequencies.Within the QHA,the Helmholtz free energy of a crystal is expressed as

whereEstatic(V,0 K) is the first-principles zero-temperature energy of a static lattice at volumeV.The second term in Eq.(1),Ezpe(V,0 K),is the zero-point energy that can be written as

The last term is the phonon free energy due to lattice vibration,given by

whereωi,q(V)is the volume-dependent vibrational frequency ofi-th mode at wavevectorq,kBand ¯hare the Boltzmann and the reduced Planck constant, respectively.Hence the vibrational entropy can be obtained by

Then the heat capacity can be calculated from vibrational entropy

The Helmholtz free energiesF(V,T) shown in Eq.(1) were calculated at 9 volume points of each specified temperature.Furthermore, the Vinet equation of state was used to fit the nine ofF(V,T) at eachTvariable within a specified temperature window considered,and the equilibrium volume at eachTwas obtained corresponding to the minimum value of fitted curve.The Gibbs free energy as a function of temperature/pressure can be written as

wherep=∂E/∂V ≈0 under the atmosphere.Transition temperature can be obtained by tracking and comparing Gibbs free energies of M and R phases as a function of temperature,satisfying Eq.(6).When Gibbs free energies of R and M phases are equal,the transition occurs naturally and theTccan be obtained simultaneously.

3.Results and discussion

3.1.Structural properties and band gap

At a temperature belowTc,the VO2is a monoclinic(M)structure with space-groupP21/c(No.14) with the unit cell parameters ofa ≈5.388 ˚A,b ≈4.614 ˚A,c ≈5.441 ˚A.The V4+ions deviates from the vertex and center of the unit cell,which results in the twist VO6octahedral cages as basic building blocks (shown in the left part of Fig.1(a)).Another feature of M phase is that the vanadium atoms title and dimerize to form zigzag-like chains with alternating different V–V bond distances, as shown in right panel of Fig.1(a).As shown in Fig.1(b),at ambient pressure and aboveTc,the thermodynamically stable phase of VO2belongs to theP42/mnm(No.136) space group with tetragonal rutile (R) structure.It has the lattice constants ofa=b ≈4.627 ˚A,c ≈2.797 ˚A.In the R-VO2phase,the V4+ions occupy the apex and body center of the tetragonal crystal cell, and six O2−ions surrounding one V4+ion constitute an octahedral VO6unit (see the left panel of Fig.1(b)).In addition,the octahedral unit forms edge-sharing chains along thecaxis.As illustrated in right panel of Fig.1(b),vanadium atoms form equally spaced V–V straight chains incdirection.

Fig.1.Crystal structure of VO2 in(a)monoclinic(M)phase with alternating V–V dimerization forming zigzag-like chains and (b) rutile (R)phase with equidistant straight V–V chains.The unit cells are indicated by black lines.Vanadium and oxygen atoms are represented by blue and red spheres,respectively.

Table 1.Band gaps and V–V distances of M and R phases by DFT+U with different U values,and compared with the experimental results.

As described above, one of the distinguishing features of VO2is its strongly correlated 3d electrons.Although the periodic plane-wave DFT calculations with typical exchangecorrelation functions, such as LDA or GGA, have significant success in describing many materials,they fail to treat strongly correlated systems.Usually,they may underestimate the band gaps and sometimes even predict qualitatively incorrect metallic ground state,e.g.,insulating nature M-VO2.This discrepancy is mainly due to inaccurate treatment of electron correlation effects.Therefore,in order to overcome this shortcoming,the DFT+Umethod with a moderate computational cost was used in the present work to better understand the electronic properties of two VO2phases.The strong Coulomb interactions between V-3d electrons are represented by the HubbardUparameter.In the case of M-VO2,the band should open a reliable gap and restore insulating character.However,previous studies have shown that the calculation results depend on the magnitude of effectiveUvalues, so we systematically examined the effects of differentUvalues on band gaps.Here,we mainly focus on the effect of HubbardUparameter on electronic structure rather than the geometrical structure.Hence,the V–V distances for the two phases shown in Table 1 are obtained from the fully optimized crystal structures by the PBE function without considering the HubbardUparameter.The PBE function can correctly reproduce two different dimerized V–V bond distances for the monoclinic phase (2.50 ˚A and 3.19 ˚A compared to the experimental values 2.62 ˚A and 3.16 ˚A,respectively).In the rutile phase,the optimized equal V–V bond length is～2.80 ˚A, which is close to the experimentally reported value of 2.85 ˚A.TakingUas the adjustable parameter in a range between 1 eV and 4 eV, the DFT+Uresults reproduce the insulating nature of M phase.The band gaps are 0.21 eV,0.44 eV,0.71 eV,and 1.00 eV withU=1 eV,2 eV, 3 eV, 4 eV, respectively.It shows a positive correlation between the band gaps and theUvalues.Meanwhile,the metallic character of R structure is preserved.From these results, we may speculate that the tuningUin the region of 2–3 eV may be required to match the optical gap (0.67 eV)from the experiments.However, one should remember that the experimental gap is measured by photoemission in an exciting progress, in which the interacting electron–hole pair is also included.A good agreement between calculation and experiment can only be achieved by considering the excitonic effects,especially in the case that there is a semiconductor or an insulator.The Bethe–Salpeter equation under the GW approximation (GW-BSE) may be appropriate for this process, but it makes the calculation expensive and cumbersome.Therefore, we would like to use a relatively inexpensive DFT+Umethod to capture the correct qualitative features of VO2system in the present work,where the reliable V–V bonds for both the phases and a finite band gap of M phase(but not the same value as the experimental one)should be reproduced.

3.2.Phonon dispersion and phase transition temperature

According to the workflow of QHA method, the phonon spectra based on force constants obtained by the densityfunctional perturbation theory (DFPT) are required prior to the actual thermodynamically physical quantities calculation.For instance, theωi,qis used to calculate zero-point energy and free energy of lattice vibration, as seen on the right-hand sides of Eqs.(2) and (3).Hence, we have first calculated the phonon dispersion curves for both M-VO2and R-VO2using DFT+UwithU= 1 eV (Fig.2).In a cell withn-atoms,there would be 3 acoustic branches, 3n −3 optical branches.There are 36 and 18 branches in the primitive cells of M and R phases,respectively.As shown in Figs.2(a)and 2(b),there is no imaginary frequency existing in the whole Brillouin zone for both the phases.It provides a direct proof of their dynamical stability under the simulated condition.From the projected phonon density of states of both the phases shown in Fig.2,the low frequency part is mainly contributed by V atoms (black solid line), and the high frequency is mainly originated from O atoms (red solid line).In addition, V atoms is more dominant for phonon modes with frequencies lower than 5 THz.These results are in agreement with the other phonon calculations given by Leeet al.[62]For the R phase, our results are different with the results from Huaet al.,[63]where there were a lot of imaginary frequencies without considering a HubbardU.This difference highlights the necessity of considering the correlation effect in rutile phase VO2.Additionally, the three non-degenerated acoustic branches, one longitude and two transverse branches, indicated that M and R phases are both anisotropic.The phonon dispersion curves by DFT+Uwith otherUvalues were also calculated.Though the details may vary slightly,there is no imaginary frequency in phonon spectra for the two phases.To keep the present text succinct,only the results ofU=1 eV is presented in Fig.2.

Fig.2.Phonon dispersion along high symmetry lines and the phonon density of states for VO2 in(a)M phase and(b)R phase.

Fig.3.Helmholtz free energy with respect to different volumes of VO2 in (a) M phase and (b) R phase at temperature window between 0 and 700 K with steps of 100 K.The filled circles denote calculated Helmholtz free energy with QHA of nine volumes at specified temperature, and the solid curves show the fitted functions according to the Vinet equation of states.The crossed points of solid and dashed line are the minimum of respective fitted functions and simultaneously represent the equilibrium volumes at each temperature.Dashed lines are a guide to the eyes.(c)The Gibbs free energies of varied temperatures for M and R phases, evolving from the dashed lines in(a) and (b), respectively.

The thermodynamic properties are calculated consequently,including the entropy,free energy,heat capacity,and so on.Figures 3(a) and 3(b) show the volume-dependent quasi-harmonic Helmholtz free energies of M and R phases.The free energies as a function of lattice cell volume at each specified temperature from 0 K to 700 K with a step of 100 K calculated with QHA forU=1.0 eV are depicted by filled circles.These calculated points are fitted by the Vinet equation of states(EOS)[64]to get the equilibrium volumes at different temperatures.As shown in Figs.3(a) and 3(b), using QHA makes it possible to preserve part of the finite temperature effect caused by the changes of phonon frequencies with the increase in equilibrium lattice volume at different temperatures.Within the framework of the QHA, the Gibbs free energy at given temperatures and a certain pressure can be obtained(the crossed points of black or red dashed line and colored solid lines).The Gibbs free energy of both M and R phases as functions of temperatures are presented in Fig.3(c).It is well known that the Gibbs free energies of the two phases at the transition point are equal.Therefore, the transition temperatureTccan be evaluated from the cross point of two Gibbs free energy curves.However, from Fig.3(c)we can find that the Gibbs free energy of R-VO2is always lower than that of M-VO2in the whole temperature window, which is counterintuitive and no obviousTc.This issue may be due to the irrationally selected HubbardU=1.0 eV,which is inadequate to correct electron–electron strongly correlated interactions.

Keep the above results in mind, one may expect to enhance theTcby adjusting theUparameter.We thus go one step further to obtain the Gibbs free energy withU=2.0 eV.As illustrated in Fig.4(a), the Gibbs free energy of M phase is lower than that of R phase in the region of 0–552.1 K,which implies that the M phase is more stable than the R phase and the phase transition occurs at 552.1 K.However theTcis higher than the experimental data.[12]Although the cross point leftwards with the decreasing HubbardUto 1.5 eV,theTc≈396.8 K [Fig.4(b)] slightly deviates from the measured value of～340 K.Additionally,the other contributions to the free energy in Eq.(1) should be considered to adjust theTcfurther.As is well known, the temperature-dependent electronic energyEelec(V,T) [Eq.(7)], usually ignored for semiconductors,should be emphasized in metal.[65]Further calculations show that the Gibbs free energy of R phase cannot be correctly describe without taking into account the electronic energy(～0.14 eV error in 0–1000 K).The reason of this issue is attributed to the large density states at Fermi level of R phase, which yields a substantial contribution of the electronic entropy.SubstitutingEelec(V,T) forEstatic(V,0 K) in Eq.(1), the modified Gibbs energies of both the phases are shown in Figs.4(c) and 4(d).The substantial contribution of the electronic entropy (Selec) decreases the free energy of R phase and hence induces the reduction ofTc.It is noticed that the phase transition temperature is 347.9 K with theUvalue of 1.5 eV(Fig.4(d)).The effective HubbardUis determined to be 1.5 eV by successfully fitting phase transition temperature.

Fig.4.The Gibbs free energy versus temperature, which is performed in a simpler and efficient framework of DFT+U combined with QHA for both M(blue dashed lines)and R phases(red solid lines)at(a)U =2.0 eV and(b)U =1.5 eV without considering electronic energy in Helmholtz free energy,(c)U =2.0 eV and(d)U =1.5 eV with electronic energy in Helmholtz free energy.The boundary of blue and red shaded region denotes the transition point from low-temperature M phase to high-temperature R phase.The corresponding temperature of boundary is the Tc,as shown in the upper right.

Fig.5.The calculated thermal properties using DFT+U combined with QHA at U =1.5 eV,including Helmholtz free energy,heat capacity, and entropy for M phase with respect to temperature, shown by orange filled circles, with red filled stars, and black solid line,respectively.The orange and red solid lines connecting filled symbols are guides for the eyes.The blue dashed line denotes the fitted heat capacity as a function of temperature by Debye approximation.Band structures and corresponding partial density of states(PDOS)calculated by DFT+U with U =1.5 eV for(b)M phase and(C)R phase of VO2.Black dashed line depicts the Fermi level.

3.3.Thermal properties and band structures with determined U

In practice, an optimal HubbardUvalue is available to correctly describe fundamental properties of both the VO2phases.In order to check the efficiency and reliability of the determinedUmentioned above, we choose to examine the thermal properties of M-VO2and the electronic band structures of both the phases with theUvalue of 1.5 eV as the examples.The calculated temperature-dependent free energy, heat capacity and entropy are obtained and presented in Fig.5(a).The heat capacityCVsimulated by the quasiharmonic approximation (the red filled stars with guided line in Fig.5(a))can be fitted according to the following equation deducted from the Debye approximation:

The parameterTDis fitted to be 728.4 K, agreeing with the experimental value of 750 K.[66]This good agreement provides a direct proof that the QHA method and theUparameter determined byTcdemonstrate the good performance for evaluating thermal properties of VO2.The electronic band structures of both the M and R phases are also presented in Figs.5(b) and 5(c) to further verify the efficiency of the determinedU.Obviously, the energy zero is shifted to the Fermi energy position in both the figures.The DFT+UwithU=1.5 eV give a～0.3 eV indirect band gap for the insulating M phase, and characterizes the metal nature of R phase.The shapes of band curves near Fermi level are comparable to the previous DFT results.[34,49]The partial density of states(PDOS) shown in the right panel of Figs.5(a) and 5(b) provide more clear evidence of the instinct characteristic of M phase and R phase.As shown in the right panel of Fig.5(b),the partial density of states of monoclinic VO2is separated near the Fermi level by a gap of～0.3 eV,which agrees with the property of its band structure discussed above.The top of valence-band sates predominantly is made up of V atoms contributions, while significant components from V atoms also contribute mainly to the lower conduction-band states.The PDOS of rutile VO2is shown in the right part of Fig.5(b).It can be noticed that Fermi level is located in the conduction bands, which indicates the metallic nature of R phase.The lower group of states below the energy−1.5 eV is mainly originated from the O atoms.The states crossing the Fermi level in the energy ranges between−0.5 eV and 2 eV are dominated by the states of vanadium atoms.It should be kept in mind what really matters is the nature of insulating M phase and metallic R phase, but not the exact value band gap.One cannot expect that the electronic band gap of the ground state always matches with the optical gap, especially in the strong correlated materials.In a word,the magnitude of the resultingUdetermined byTcbased on QHA is physically reliable.

4.Conclusions

Using first-principles calculations with the quasiharmonic approximation, the temperature-dependent Gibbs free energy curves of insulating M phase and metallic R phase of VO2with a certain HubbardUcan be plotted.The transition temperatureTcis determined by identifying the cross point of Gibbs free energy curves during the phase transition.By comparing the calculatedTcwith the experimental reference value, the effective HubbardUvalue is turned out to be 1.5 eV for VO2.Moreover,the calculations with thisUvalue may correctly capture the metallic and insulating properties of the M and R phases, respectively.In conclusion, this work demonstrates that the first-principles calculations in combination with QHA can be used to determine the effective HubbardUvalue of VO2by analyzing the phase transition temperatureTc.This approach is promising for determiningUof other strongly correlated materials.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant Nos.21933006 and 21773124)and the Fundamental Research Funds for the Central Universities Nankai University(Grant Nos.010-63233001,63221346,63213042, and ZB22000103).K.L.acknowledges the support from the China Postdoctoral Science Foundation (Grant No.2021M691674) and the Hefei National Laboratory for Physical Sciences at the Microscale(Grant No.KF2020105).