GGA-1/2 self-energy correction for accurate band structure calculations: the case of resistive switching oxides

First-principles calculation has become an indispensable methodology in revealing the working principles of nanoscale electronic devices, but ultra-large supercells are usually required in modeling the devices with critical metal/dielectric interfaces. Traditional density functional theory within the generalized gradient approximation (GGA) suffers from the inaccurate band gap problem when metal oxides are present, but they serve as the core component in resistive random access memory (RRAM), which is a promising path for novel high speed non-volatile memories. To obtain improved oxide band gaps, we applied the efficient GGA-1/2 method for self-energy correction, whose computational load is at the same level as standard GGA. In particular, we have investigated the influence of exchange-correlation functional flavors on the GGA-1/2 band structures, taking four important binary oxide RRAM materials (α-Al2O3, r-TiO2, m-ZrO2 and m-HfO2) as benchmark examples. Five GGA functionals (PBE, PBEsol, PW91, revPBE and AM05) were considered and their band structures were compared in detail. We have found that the performance of GGA-1/2 is comparable to state-of-the-art GW and generally superior to the HSE06 hybrid functional. Among the five GGA functionals, PBEsol yields the best results in general. In addition, the applicability of a single self-energy potential for various GGA-1/2 flavors is discussed. Our work provides a guide to the GGA flavor selection, when applying the GGA-1/2 method to metal oxides.


Introduction
With the rising of big data [1], cloud computing [2], internet of things [3] and other technologies in recent years, an exploding amount of information requires storage and analysis, for which fast, low power consumption, high density and non-volatile memories are suitable. However, the NAND FLASH memory that dominates the present non-volatile memory market suffers from intrinsically slow erasing operation, opening an urgent demand for next-generation non-volatile memories [4,5]. Resistive random access memory (RRAM), which achieves data storage through reversible resistance changes in certain dielectrics, is one of the leading candidates for this purpose. Reproducible resistive switching phenomena have been widely observed in oxides [6][7][8][9], chalcogenides [10] and even organic materials [11]. Due to their simple device structure (metal/insulator/metal), high compatibility to standard CMOS flow line and the prospect of 3D integration, certain binary oxides such as Al 2 O 3 [12], HfO 2 [13,14], TiO2 [15,16], and ZrO 2 [17] have attracted more attention. In recent few years, HfO2 has been regarded in some works [18] as the most promising candidate material for commercial application.
Although oxide-based RRAM has experienced a rapid development for more than a decade, it has not reached large-scale industrial applications yet. There are still some fundamental issues to be resolved in terms of the resistive switching mechanism, reliability, integration schemes, and so forth [15]. For instance, the mesoscopic physical mechanisms of switching and the corresponding models, without which a proper device optimization cannot be achieved, still require more in-depth study. Most binary metal oxides are believed to manifest a transition from a high resistance state (HRS) to a low resistance state (LRS) upon the formation of conductive filaments (CFs) connecting the two electrodes. Yet, at present there are still controversies on the exact composition and geometry of the CFs [19][20][21][22]. Since it is difficult to directly observe and determine experimentally the exact composition of the CFs within the dielectric media, theoretical approaches can be particularly helpful.
Among atomistic simulation techniques, ab initio calculations based on density functional theory (DFT) [23] are highly successful in the prediction of total energies, bond lengths, and vibration frequencies. However, DFT within local density approximation (LDA) or generalized gradient approximation (GGA) severely underestimates the band gaps of semiconductors and insulators, which limits its capability of predicting the band alignment in metal/semiconductor contacts, a significant task in understanding charge transmission through CFs [24]. Various theoretical approaches have been proposed to overcome this shortcoming, thus providing more reliable predictions of physical properties that depend on excited states. Hedin's GW approximation is considered to be the state-of-the-art, which is still the most accurate method for electronic structure calculation. It calculates the quasiparticle energies in terms of the perturbation theory [25], and is a beyond-DFT technique. Hybrid functionals [26], which combine standard DFT and Hartree-Fock exchange functionals through a mixing parameter, are also well-known and robust alternatives, though they sometimes require an adjustment of the mixing parameter. Both hybrid functional and GW methods demand much greater computational efforts than standard DFT-LDA or DFT-GGA, setting severe limitations to the size of the simulated systems. Indeed, to simulate realistic RRAM devices, large supercells with hundreds of atoms are typically required. For instance, to investigate the atomic structure of a CF through the dielectric, a large crosssection area is indispensable to avoid undesired interactions between a CF and its periodic images [13]. Therefore, realistic ab initio RRAM device calculations require both a high computational speed and relatively accurate band gaps. Traditionally, one resolves the band gap problem in DFT-LDA/DFT-GGA by employing LDA+U and GGA+U [27][28][29][30], or resorts to hybrid functionals [31]. While the adoption of hybrid functional increases the computational load considerably, LDA+U is almost as fast as conventional LDA. However, LDA+U is not suitable for many of these RRAM oxides. For example, Al 2 O 3 , TiO 2 , ZrO 2 and HfO 2 are all ordinary oxides without strong electron correlation, in which case LDA+U does not work well unless an additional large U p correction to oxygen is considered [27], reflecting the fact that the band gap inaccuracy lies primarily with the oxygen anion rather than the metal cation.
The LDA-1/2 method corrects the spurious electron self-interaction term in LDA using the Slater halfoccupation technique [32]. As the conduction band electron usually occupies Bloch-like states with nearly-null self-energy, corrections are only carried out to the holes that are typically localized on the anions. Instead of calculating explicitly the self-energy of the hole, such correction can be achieved by simply introducing a socalled self-energy potential in real space, which acts as an external potential to the region where the hole resides. The exact form of self-energy potential for any atom can be derived ab initially from atomic calculations using DFT, with the desired exchange-correlation (XC) flavor. Such self-energy potential for the anion is then trimmed to a finite range to avoid overlapping between neighboring anions. Going one step further from LDA-1/2, a natural extension is GGA-1/2, which was previously shown to bring about improvement over LDA-1/2 in some cases [33]. While LDA tends to underestimate the lattice constants, in the framework of GGA there are certain XC functionals like PBEsol [34], Wu-Cohen [35] and AM05 [36] that can predict rather accurate lattice constants. Hence, GGA-1/2, which corrects the electronic band structure calculated with GGA, ought to be a better computational method for RRAM simulation than LDA-1/2. However, to our best knowledge there is no theoretical study on the influence of GGA flavors towards the accuracy of GGA-1/2 band structures. Hence, from both theoretical and application aspects, it is of great significance to carry out a systematic study on the GGA-1/2 band structures built upon various GGA flavors.
In this work, we adopt four relevant binary oxide RRAM materials as benchmarks. Five GGA XCs are considered and their band structures are compared in detail. Our results offer guidance in properly choosing GGA flavors to achieve efficient and accurate GGA-1/2 band structure calculations for the four binary oxide semiconductors, i.e. corundum Al 2 O 3 , rutile TiO 2 , monoclinic ZrO 2 and monoclinic HfO 2 . Our results can be readily extended to more oxides.

Crystal structures
We studied the room temperature phases of Al 2 O 3 , TiO 2 , ZrO 2 and HfO 2 , namely corundum (α-Al 2 O 3 ), rutile-TiO 2 (r-TiO 2 ), monoclinic ZrO 2 (m-ZrO 2 ) and monoclinic HfO 2 (m-HfO 2 ), shown in figure 1. Corundum α-Al 2 O 3 crystallizes with trigonal symmetry in the space group R c 3 , where oxygen atoms are arranged by a slightly distorted hexagonal close packing, in which two-thirds of the gaps between the octahedra are occupied by aluminum ( figure 1(a)) [37]. Rutile TiO 2 presents a tetragonal unit cell with the titanium atoms surrounded by an octahedron of six oxygen atoms. The oxygen coordination number is 3, resulting in a trigonal planar geometry, as illustrated in figure 1(b) [38]. The ground state monoclinic phases of ZrO 2 and HfO 2 are isomorphic, both with the baddeleyite structure (P2 1 /c) as shown in figure 1(c) [39]. Unlike TiO 2 , which features six-coordinated Ti in all phases, m-ZrO 2 and m-HfO 2 consist of seven-coordinated zirconium or hafnium centers. This discrepancy is attributed to the larger size of Zr/Hf atoms relative to the Ti atom.
2.3. GGA-1/2 for improved band structures GGA-1/2 generalizes to solids the original idea of Slater's half occupation scheme, which was initially designed for ionization energy calculation in isolated atoms. Its formalism starts from Janak's formula, which offers an interpretation of the DFT single particle eigenvalues e α as the variation of the total energy E with respect to its Assuming a linear relation between the eigenvalue and its occupation [48], it is straightforward to show that Equation (2) means that the energy difference between the ground state energy E(0) and the total energy of the ion E(−1) with one electron removed from state α, is given by the energy eigenvalue of state α and a so-called "self-energy term" S α . The analytical expression for S α [32], dr dr n r n r r r higher order terms 3 3 3 indeed has the form of a self-energy, which renders its name. If one regards S α as the quantum mechanical average of a "self-energy potential" V s , it follows that V s can be obtained as the difference from the DFT-calculated electrostatic potentials between the neutral atom and the half-ionized ions: As mentioned in section 1, the self-energy potential must be trimmed by some function Θ(r) [32]. To obtain the rectified band gaps for solid state materials, V s is added to the pseudopotentials of the anions and a second self-consistent run using the modified pseudopotentials is carried out upon the same relaxed geometry. Since the self-energy potentials serve as external potentials in the solid, the ground state energies given by GGA-1/2 selfconsistent calculations are not physically meaningful themselves. Nevertheless, band diagrams and density of states are properly recovered by GGA-1/2 [32]. Hence, in our work all total energies were taken from GGA calculations, while the band diagrams were taken either from GGA (uncorrected) or from GGA-1/2 (corrected) calculations.
Since all the materials under investigation are binary metal oxides, the self-energy correction ought to be done only for oxygen anions. The optimal cutoff radius in Θ(r) to trim the corresponding V s (with (1/2)e subtracted from the O 2p orbital [32]) was obtained variationally upon maximizing the band gap [32]. Electronic structures were obtained with this optimal cutoff radius only, and no empirical parameters were involved.
The self-energy potentials with all these GGA XC flavors (PBE, PBEsol, PW91, revPBE, AM05) were derived from atomic calculations using a modified ATOM code (supplied with the Siesta simulation package [49]). These self-energy potentials do not depend on the ab initio simulation code, and can be regarded as only specific to the XC flavor. We have also created a web-based self-energy correction program where pseudopotentials modified through the inclusion of GGA-1/2 self-energy potentials for all these XC flavors can be generated online, at http://www.eedevice.com/dft-half.

Geometric structures and partial density of states
The optimized lattice parameters for these four crystals are summarized in table 1, which also shows the corresponding experimental values [50][51][52]. Due to their intrinsic similarity, PBE and PW91 yield very similar lattice parameters, while the lattice parameters obtained with PBEsol and AM05 also share a similar trend. The lattice parameters obtained with revPBE are always the largest among all considered XC functionals. Compared with experimental values, PBEsol and AM05 are closest to data, while PBE and PW91 usually predict slightly larger lattice constants. On the other hand, revPBE yields the worst agreement of all functionals tested.
The partial density of states for α-Al 2 O 3 , r-TiO 2 , m-ZrO 2 and m-HfO 2 obtained with GGA-PBE are shown in figure 2. In any material the valence band maximum (VBM) mainly consists of O-2p states, which confirms that the hole is located around oxygen anions and the self-energy correction should indeed be done to oxygen only.

Band structure results based on optimized lattice parameters
The band gaps calculated with PBE, PBEsol, PW91, revPBE, AM05, HSE06 and GGA-1/2 (including the five different XC flavors) for the four compounds are listed in table 2, where fully optimized lattice parameters were utilized for each XC. As expected, DFT-GGA severely underestimates the energy gap. All band gaps obtained with conventional GGA are far below the experimental values by 1.16 eV∼2.67 eV. Among the four materials, α-Al 2 O 3 is a wide band gap insulator with the VBM and conduction band minimum (CBM) both located at Γ, which is consistent with previous results [53,54]. The largest obtained band gap for α-Al 2 O 3 based on DFT-GGA is 6.13 eV, using the AM05 XC. However, that is still 2.67 eV smaller than the typical experimental value (8.8 eV) [53], and even smaller than the previous GW results (9.36 eV and 9.78 eV) [55,56]. After self-energy correction, the GGA-1/2 band gaps obtained with the five XC functionals are all above 8 eV (table 2), which is comparable to the HSE06 values 8.09 eV (calculated previously for a PBEsol-optimized structure using normconserving pseudopotentials [54]) and 7.70 eV (calculated in this work for a PBE-optimized structure).    [57][58][59][60][61]. Compared with experimental band gap values [58][59][60][61], the GGA band gaps differ by about 1.16 eV∼1.56 eV (table 2). Yet, the GGA-1/2 band gaps are quite satisfactory, agreeing well with experimental values by showing merely −0.1 eV to 0.4 eV mismatches. In addition, our GGA-1/2 results are also very close to the HSE06 (3.32 eV and 3.39 eV [58]) and GW (3.46 eV) [58] results.
Regardless of the computational method, both m-ZrO 2 and m-HfO 2 are predicted to be indirect band gap semiconductors with the VBM and CBM located at Γ and B, respectively. In general, the calculated GGA band gaps are lower than typical experimental values by about 1.30 eV∼2.32 eV (for m-ZrO 2 ) and 1.80 eV∼1.94 eV (for m-HfO 2 ). After self-energy correction, the GGA-1/2 band gaps differ from the experimental value by merely −0.59 eV∼0.47 eV for m-ZrO 2 and −0.27 eV∼−0.06 eV for m-HfO 2 , respectively. Regarding higher order calculations, the HSE06 band gap for m-ZrO 2 (5.20 eV in this work and 5.14 eV from [63]) is lower than our GGA-1/2 results as well as the available experimental values. Moreover, our GGA-1/2 results are comparable to that from GW 0 (5.34 eV [62] and 5.42 eV [64]) and larger than that of G 0 W 0 (4.99 eV) calculations [62]. Comparing the five XC flavors, PBEsol yields the best band gap that differs from experimental by less than  Figure 3 shows the best (among all XC functionals) GGA-1/2 band structures for each material, along with conventional DFT-GGA and HSE06 results. The GGA-1/2 band structures are in line with DFT-GGA and HSE06 in terms of the morphology and the intrinsic excitation characteristics (such as the locations of VBM and CBM). This confirms that for these oxides, the good accuracy of GGA-1/2 is not limited to the fundamental gap, but also for the overall band structures, where GGA-1/2 is comparable with computationally intensive methods such as hybrid functionals.

Band structure results based on experimental lattice parameters
Since each GGA XC leads to a specific set of optimum lattice parameters, the band diagrams calculated in section 3.2 were not based on a unique unit cell for each material. To better evaluate the difference among various XCs, we further calculated the band gaps of α-Al 2 O 3 , r-TiO 2 , m-ZrO 2 and m-HfO 2 with both GGA and GGA-1/2 using their experimental lattice parameters. Yet, to avoid internal forces, the atomic coordinates within the unit cell were fully relaxed. The calculated band gaps for all five XC flavors are listed in table 3, and figure 4(a) compares those band gaps obtained with optimized lattice parameters with fixed experimental lattice parameters. Among the five XCs, the obtained band gap of r-TiO 2 with experimental lattice parameters is slightly smaller than that of optimized structure for PBEsol, PW91 and AM05, regardless of using GGA or GGA-1/2. On the other hand, the results of PBE (both GGA and GGA-1/2) and revPBE (only GGA-1/2) are slightly larger than that of optimized structure. In addition, for α-Al 2 O 3 the band gaps are much improved by fixing to experimental lattice constants, especially for the revPBE XC whose GGA-1/2 band gap now only deviates from the experimental value by 0.15 eV. It is therefore clear that the original poor GGA-1/2 band gap using equilibrium lattice constants is strongly related to the over-estimation of lattice constants. For instance, the equilibrium c-axis lattice constant was predicted to be 13.211 Å using the revPBE XC, which is 1.7% larger than the experimental value. Such improvement also implies that using correct lattice constants is more important than minimizing the stress, in order to obtain band gaps close to measured values. Finally, we notice that the lattice constant effect on m-ZrO 2 and m-HfO 2 is relatively tiny, because the band gaps derived from either optimized lattice or experimental lattice differ by no more than 0.1 eV.
It is well-known that AM05 and PBEsol yield the most accurate lattice parameters among common GGA XCs. Therefore, it is natural to infer that for these two XCs the GGA-1/2 band gaps do not vary much when switching to experimental lattice parameters. However, the quantitative relation is still unclear yet. Hence, we have also listed explicitly in figure 4 lattice)] using GGA-1/2. Notice that while the differences for m-ZrO 2 and m-HfO 2 are slight regardless of the XC used, for α-Al 2 O 3 and r-TiO 2 the differences can be a lot more pronounced. Moreover, among the five GGA XCs the ΔE g values obtained with PBEsol have the least fluctuation in general, slightly better than AM05 (both PBEsol and AM05 are considerably better than the other three GGA XCs).

Influence of different pseudopotential forms
Pseudopotentials of different forms can have great influence on the calculated electronic structures, especially for transition metals due to the complexity of d electrons. The pseudopotential of Al is sufficiently simple because only 3 electrons need to be considered as in the valence. On the other hand, for the 4d metal Zr, the pseudopotential of the Zr_sv form (including 4s, 4p, 4d and 5s electrons in the valence) is unavoidable in VASP. There is thus little room to test the influence of pseudopotential forms for Al and Zr. Hence, we only considered different GGA-PBEsol pseudopotential forms for Ti and Hf available with the VASP code, which are Ti(Hf) (four valence electrons, 4e for short), Ti(Hf)_pv (ten valence electrons, 10e for short), and Ti(Hf)_sv (twelve valence electrons, 12e for short). As shown in table 4, the calculated band gap of r-TiO 2 is much more sensitive to the choice of pseudopotential than that of m-HfO 2 . For r-TiO 2 , at least the Ti_pv (10e) pseudopotential should be used in order to get accurate results. In addition, all the three pseudopotential forms considered here for Hf are suitable for various calculations. Thus, on account of both computational accuracy and efficiency, in carrying out GGA-1/2 electronic structure calculations we recommend Ti_sv (12e) or at least Ti_pv (10e) pseudopotentials for r-TiO 2 , while for m-HfO 2 all three available pseudopotential forms are suitable.

Universality of the GGA-1/2 self-energy potentials
In all the GGA-1/2 calculations above, we have derived the self-energy potentials by subtracting the potential of a charged ion with (1/2)e missing, from that of a neutral atom. Both neutral and ionized atomic potentials were calculated ab initially using DFT-GGA for each XC flavor. These self-energy potentials in fact do not differ much since the subtraction may eliminate most of the specific features belonging to each specific XC flavor. To investigate whether a unique self-energy potential can be used in the GGA-1/2 calculation with different GGA XCs, we have manually listed PBE and PBEsol as the baseline XCs. This means that we subsequently carried out GGA-1/2 calculations with all five XCs, but utilizing the set of self-energy potentials only from the baseline XCs. For all oxides, the corresponding band gap results are shown in table 5. Compared with table 2, the maximum GGA-1/2 band gap variation upon switching to the PBE-derived oxygen self-energy potential, is merely 0.02 eV among all XCs and all oxides. When exclusively using PBEsol-derived oxygen self-energy potential, the maximum band gap variation is enlarged to 0.06 eV, which is still a tiny value. We conclude that for these materials the self-energy correction is not significantly dependent on the XC flavor, and the PBE self-energy potentials may fit other XC flavors in carrying out GGA-1/2 calculations.
3.6. Influence of the self-energy potential trimming scheme In the GGA-1/2 calculations, a cutoff function Θ is employed to trim the self-energy potential: where p is a power index that should be sufficiently large to ensure a sharp trim, and r cut is the cutoff radius. In the original work of Ferreira et al the value p=8 was chosen to avoid possible convergence issues, but higher values can also be used. Indeed, we re-calculated the band gaps for the four oxides using GGA-1/2 with p=20, and encountered no issue related to numerical convergence. Compared with p=8, the GGA-1/2 band gaps with p=20 (listed in table 6) increase by about 0.2 eV (for α-Al 2 O 3 ) and 0.09 eV∼0.14 eV (for r-TiO 2 , m-ZrO 2 and m-HfO 2 ). The mean differences between the calculated band gap values and the experimental results are diminished as p increases from 8 to 20, indicating that p=20 is more favorable for these calculations because it includes the self-energy more accurately. In order to explain this phenomenon, in figure 5 we plot the real-space spatial distribution of the hole (residing at the VBM in k-space) and the electron (residing at the CBM in k-space) for r-TiO 2 as an example. Since the hole and electron overlap, when the p value is low the self-energy potential reaches the electron density, lowering the CBM as well as the VBM. On the other hand, as p increases from 8 to 20, the self-energy potential acquires a sharper boundary, resulting in less overlap with the electron distribution, which yields larger band gaps. In the meantime, we notice a variation of the optimal oxygen self-energy cutoff radius when using different p values. As plotted in figure 6, the cutoff radii are in general reduced by around 0.2 a.u.∼0.3 a.u. when p is increased from 8 to 20. Such phenomenon is common for all the materials under investigation. However, as long as p is fixed, the optimal r cut is relatively constant across these oxides, showing good transferability of the self-energy potentials.  The hybrid functional HSE06, depicted as hollow triangles in figure 7, clearly shows great improvement over conventional GGA, but its band gaps are generally smaller than those from GW calculations. On the other hand, the GGA-1/2 band gaps are generally better than HSE06 results, except for r-TiO 2 . Among the five XC flavors, PBEsol demonstrates the best performance in general, in accordance with our previous results discussed in Section 3.3.

GGA-1/2 in the study of RRAM devices
After demonstrating the advantages of GGA-1/2 and identifying PBEsol as the XC with the best accuracy, we here list some possible applications of the GGA-1/2 method in investigating some fundamental problems in RRAM, especially regarding the nature of the CF: • Prediction of CF composition and structure i n binary oxide-based RRAM cells. It is commonly agreed that oxygen vacancies (V o 's) are key to CF formation. However, for each material it is far from certain what amount of V o should render CF formation with the Table 6. Calculated GGA-1/2 band gaps (unit: eV) for α-Al 2 O 3 , r-TiO 2 , m-ZrO 2 and m-HfO 2 with different exchange-correlation functionals, with a power index p=20 in the self-energy potential trimming function. The differences between band gaps obtained with p=8 and p=20 (ΔE g ) are also listed.    typical resistive switching behaviors found in experiments [70]. On the one hand, in the study of crystallinephase sub-oxides [7], their conductive nature in principle ought to be predicted only by higher level techniques such as hybrid functionals and GW [24], since LDA/GGA incorrectly predicts the band gap, which tends to produce incorrect I-V characteristics. LDA/GGA can even predict a semiconductor to be a metal as it does to Ge [71,72]. On the other hand, to evaluate more complex structures, such as amorphous phase CF candidates with various V o concentrations, or CFs located in grain boundaries [73,74], hybrid functionals and GW are computationally too demanding to be used. As we have shown, GGA-1/2 can efficiently deal with bulk crystalline oxides, as well as their amorphous phases or various interfaces.

PBE
• Investigation of the horizontal scaling limit of RRAM cells . The horizontal scaling limit means the minimum capacitor area that can sustain reliable LRS and HRS for a given RRAM technology. To this end one needs to establish the minimum stable area of a CF, and the minimum area of the insulating oxide phase surrounding the CF. The capability of maintaining the conductive state is important for the CF in small area capacitors, while the thermodynamic stability of the CF is another key factor. Hence, a CF-inside-oxide composite supercell should be studied using ab initio techniques. Initially we have used GGA to study the size limit of Hf suboxides inside HfO 2 , even though the band gap was inaccurate [13]. Recently we have also investigated the area limit of Hf CFs inside HfO 2 , where GGA-1/2 was implemented such that the band gap of HfO 2 was correctly recovered [75].
• Investigation of the vertical scaling limit of RRAM cells . In the vertical direction the RRAM cells also possess minimum thickness limit for the dielectric, since an even thinner dielectric may suffer from electron tunneling such that the HRS becomes leaky, which goes against the goal of a non-volatile memory. This requires setting up a 5 nm thick dielectric layer typically, encapsulated by two electrodes. Such setup can contain more than 500 atoms. Moreover, a series of dielectric thicknesses should be considered to extract the minimum dielectric thickness from the electronic structures of these capacitors. This computational task, which requires both band gap accuracy and numerical efficiency, can be addressed by GGA-1/2.
• Understanding device performance and reliability.
GGA-1/2 can also be used to investigate the reasons why a specific RRAM shows certain performance limitations. A typical example is the resistance window of the TiO 2 and HfO 2 RRAMs, where the latter is substantially larger [13,76]. Wu et al [77] calculated the Schottky barriers of Ti 4 O 7 /TiO 2 using GGA-1/2, where Ti 4 O 7 is the typical CF composition in TiO 2 -based RRAMs [7]. They showed that the contact is Ohmic, thus the HRS of TiO 2 RRAMs does not possess the necessary barrier for electron injection from the CF to the dielectric. On the other hand, the CF/HfO 2 Schottky barriers are sufficiently large, no matter whether pure Hf or conductive sub-oxides were used as the CF model. GGA-1/2 is very useful for calculating the electronic structures of these models since previous calculations either suffered from band gap inaccuracy [78], or were forced to adopt short models for the dielectric to enable heavy hybrid functional calculations [31].

Conclusions
In summary, we performed a thorough investigation into the accuracy of the GGA-1/2 electronic structure calculation method, by comparing five different exchange-correlation functionals (PBE, PBEsol, PW91, revPBE and AM05). Four metal oxide crystals (α-Al 2 O 3 , r-TiO 2 , m-ZrO 2 and m-HfO 2 ) were selected for benchmark calculations, based on their relevance in resistive memory devices, for which the simulation requires exceptionally large supercells in order to clarify the device physics. Such tasks ask for fast and accurate schemes for the calculation of band gaps and band offsets. We compared the accuracy of each GGA functional flavor and have drawn the following conclusions: (1) GGA-1/2 within any of these exchange-correlation functional flavors predicts reasonable band gaps, correcting the band gap underestimation of conventional GGA.
(2) For the oxides investigated, the GGA-1/2 band gaps fall between the HSE06 and GW results, while the GGA-1/2 band structure morphology is in agreement with the computationally-demanding GW method.
(3) Since PBEsol and AM05 predict the most accurate lattice parameters, GGA-1/2 based on these two flavors yields band gaps that are relatively insensitive to the choice of experimental or optimized lattice parameters. In particular, PBEsol performs slightly better than AM05.
(4) While the atomic potentials obtained with different GGA flavors are fairly different, the corresponding selfenergy potentials are relatively insensitive to the XC flavor that was used in the self-energy potential generation. In case a self-energy potential with a specified GGA XC is missing, one may resort to the PBE self-energy potential.
(5) In the implementation of GGA-1/2, increasing the power index of the self-energy cutoff function from 8 to 20 can further improve the calculated band gaps for typical metal oxides.
Our work shows that GGA-1/2 is an efficient and accurate method for calculating the electronic band structures for various resistive switching oxides, for all the GGA exchange-correlation functionals investigated. Nevertheless, the combination GGA(PBEsol)-1/2 stands out as the optimum scheme that offers the most accurate band gaps for these oxides.