{'fraction_non_alphanumeric': 0.03653561500415178, 'fraction_numerical': 0.019314776706740314, 'mean_word_length': 4.315678372673191, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
QSGW calculation of the work functions of Al(111) and Al(100) surfaces

Modifications to the quasiparticle self-consistent GW (QSGW) method needed to correctly describe metal/vacuum interfaces and other systems having extended regions with small electron density are identified and implemented. The method's accuracy is investigated by calculating work functions for the Al(111), Al(100), and Al(110) surfaces. We find that the results for work function do not depend on the DFT functional employed to calculate the starting Hamiltonian and that QSGW yield results in quantitative agreement with data from ultrahigh vacuum experiments.


I. INTRODUCTION
The work function is not only a most important quantity characterizing the surface of a metal, it also directly affect surface phenomena like growth rate, the form of crystallites, sintering, catalytic behavior, adsorption, chemical reactions, surface segregation, and formation of grain boundaries. In addition, the work function largely determines rates of electron surface emission and is therefore of applied interest for optimizing thermionic emitters 1 , where a low work function is sought, and pulsedpower components, e.g, in transmission lines, where a high work function is required. An ability to make accurate theoretical calculations of work functions is therefore of great interest and importance for the development of new high-performing materials.
Presently, most calculations of work functions are made within density functional theory (DFT) 2,3 using either the local density approximation (LDA) or the generalized gradient approximation (GGA). Unfortunately, the accuracy of the DFT-based methods in calculation of the work function is not always satisfactory and results depend on the DFT functional used (see, e.g., recent review of available theoretical and experimental values of work function for a number of metals in Ref. 4). For instance, the work functions of Al surfaces calculated by using the GGA with Perdew-Burke-Ernzerhof (GGA/BPE) functional 5 are approximately 0.2 eV lower then corresponding LDA values calculated with Ceperley-Alder (LDA/CA) functional 6 and approximately 0.3 eV lower then corresponding LDA values calculated with Barth-Hedin (LDA/BH) functional 7 (see Table I). This example highlights the need for methods capable of reaching 0.1 eV (or better) accuracy for theoretical prediction of the work function of metals.
The GW approximation of Hedin 8 is a well established method which yields highly accurate quasiparticle (QP) energies for bulk materials 9,10 . GW calculations are usually performed in a non-self-consistent manner by using the LDA Green's function, G, and screened Coulomb interaction, W (so-called G 0 W 0 method). The results obtained by the G 0 W 0 method depend on the quality of underlying DFT wave-functions and eigenvalues and the agreement with experimental data worsen if the DFT description is not sufficiently accurate. For example, G 0 W 0 fails to describe the band gap of NiO 11 . Recently, Faleev, van Schilfgaarde and Kotani 10,11 developed the so-called Quasiparticle Self-consistent GW (QSGW) method which is independent of DFT and demonstrated that results for QP energy levels of bulk materials obtained by the QSGW are in better agreement with experiment then results obtained by the standard G 0 W 0 method. In contrast to the G 0 W 0 method, QSGW describes correctly also strongly correlated materials like NiO or MnO 11 .
Although calculations of surface QP energies were performed already more then two decades ago 12 , GW calculations for surfaces and other non-bulk systems remains rare due to the demanding computational requirements. Most commonly, the GW method has been applied to study an image potential and corresponding image states of insulating 13 and metallic [14][15][16] surfaces and clusters 17 . Recently, the GW method was used to investigate the image potential-induced renormalization of the molecular electronic levels for molecules adsorbed on metal surfaces 18,19 and thin insulator films 20 . Note that the image potential cannot be obtained by the LDA or GGA approaches since they do not include non-local polarization effects that are present in GW theory through the W operator. Another application of the GW method to non-bulk materials that is becoming an active area of research is QP calculations of transport properties of nanoscale systems [21][22][23] .
To the best of our knowledge, two GW studies of the work functions of metals have been published to date. Morris et al 24 calculated the work functions for Al(111), Al(100), and Al(110) surfaces using the G 0 W 0 method and making a jellium approximation. They also included vertex corrections to the self-energy and W (evaluated on a homogeneous electron gas level) that resulted in a significant, over 1 eV, underestimation of calculated work functions, attributed to inherent self-interaction error. Heinrichsmeier et al 25 proposed a new non-local parametrization of the exchange-correlation functional derived from the G 0 W 0 calculations for jellium surfaces (the V xc (GW) method), and applied the method to real (111) and (100) surfaces of Al and Pt. A conclusion of both studies was that the values of the Al(111) and Al(100) work functions obtained by the G 0 W 0 method are significantly worse then corresponding LDA values as compared to the experimental data (see Table I). Since both these GW calculations involve the jellium approximation, the question remains how well the fully atomistic GW method could describe the work functions of metals as compared to the DFT results and experiment.
In this paper, we present results of work function calculations for the Al(111), Al(100), and Al(110) surfaces evaluated within fully atomistic G 0 W 0 and QSGW approaches. As we show below, the values of the work function obtained by these two methods only insignificantly deviate from each other (within 0.02 eV), do not depend on the initial DFT functional, and are in excellent agreement with experimental data. On the other hand, the DFT approaches predict different results for the work function depending on the functional used. The paper is organized into two main sections describing the methods and the results followed by a concise summary.

II. METHOD AND COMPUTATIONAL DETAILS
Let us first briefly describe the computational approach. QSGW is a method to determine nonlocal (but static and Hermitian) optimum one-particle Hamiltonian H 0 in a self-consistent way 10,11,26 . First, starting with a trial Hamiltonian H 0 (usually, the LDA Hamiltonian is used as the first-iteration H 0 ) the self-energy Σ(ω) is calculated in the GW approximation. The static self-energy is defined in the basis of the eigenfunctions ψ kn (r) of the Hamiltonian H 0 as follows 10,11,26 Σ k nn ′ = Re ψ kn |[Σ(ε kn ) + Σ(ε kn ′ )]/2|ψ kn ′ , (1) where k is the wave vector, n is the band index, ε kn denote eigenvalues of H 0 , and Re means to take the Hermitian part. Next, the Hamiltonian H 0 is updated for each iteration using Σ k nn ′ instead of the usual LDA exchangecorrelation potential where V ext is the external (nuclei) potential and V H is the Hartree potential. This procedure is iterated until self-consistency is reached. Here, self-consistency is defined as when Σ k nn ′ generated by H 0 is identical (within a small tolerance) to the Σ k nn ′ that enters into H 0 . It has been shown 10 that this procedure in an approximate way minimizes the difference between the full nonlocal, non-static and non-Hermitian GW Hamiltonian H(ω) = − ∇ 2 2m + V ext + V H + Σ(ω) and Hamiltonian H 0 (that is why it is called 'optimum'). Note that H(ω) is a functional of H 0 because both V H and Σ(ω) calculated in the GW approximation depend on eigenfunctions generated by H 0 . Hence, the iteration procedure described above self-consistently determines both H(ω) and the corresponding optimum H 0 . The G 0 W 0 method is simply the first iteration of described above cycle, selfenergy is obtained from Eq. (1) using LDA wave function and energies, and first-iteration G 0 W 0 Hamiltonian is constructed by Eq. (2). In the present work, we used the experimental lattice constant of Al at zero temperature, 4.025Å 27 . The Al(111), Al(100), and Al(110) surfaces were modeled by a (1 × 1) surface unit cell in xy-directions and a periodic combination of N A Al layers and N V vacuum layers in z-direction. The vacuum layer was of the same width as the Al layer, and contains so-called floating basis orbitals 28 placed instead of the atomic muffin-tin orbitals. N A ranged from 4 to 12 and N V ranged from 6 to 10 were used to analyze the convergence of the results with respect to these parameters. The work function, Φ, is defined as the difference between the electrostatic potential at a point far from the surface and the Fermi energy; Φ = V es (∞) − ε F . In our calculations V es (∞) was estimated as electrostatic potential in the middle of the vacuum slab.
It is known that GW calculations of band gaps in semiconductor thin films 29,30 and molecular chains 31 performed in repeated-cell geometries converge slowly with vacuum thickness due to the long-range nature of the non-local screened Coulomb interaction. Rozzi et al. 31 investigated this problem and developed a Coulomb cut-off scheme to eliminate the long-ranged slab-slab interaction. They found that the introduced Coulomb cutoff parameter mostly affects delocalized unoccupied states. On the other hand, both quantities that enter the expression for the work function: the Fermi energy and electrostatic potential are affected only by occupied states that are localized inside the metal slab and therefore only weekly depend on the thickness of the vacuum slab. Consequently, the work function converges quickly with the number of vacuum layers: we found that N V =6 is sufficient to determine the work function to 0.005 eV accuracy.
The QSGW method is implemented as an extension of the all-electron full-potential linear muffin-tin orbital (LMTO) program suite. The diagram of the selfenergy and density/potential self-consistency cycles that includes the LMTO and GW parts of the code are shown on the Figure 2 of Ref. [26]. The description of the basis sets and other details of the LMTO and QSGW implementations can be found in Refs. [26] and [28].
The surface Brillouin zone (BZ) integration in the LMTO part of the code were performed with (22 × 22) and (24 × 24) Monkhorst-Pack meshes 33 (k LMT O -mesh). The GW self-energy was calculated with (6 × 6), (8 × 8), and (10×10) meshes in the surface BZ (k GW -mesh). The modified offset-Γ method designed to treat anisotropic systems was employed to perform k-integration of surface BZ in the GW part of the code (the method is described in Ref. [26], following the Eq. (53)). The GW part of the code, where the self-energy is calculated given the eigenfunctions and eigenvalues of the H 0 generated by the LMTO part, is significantly more computationally demanding than the LMTO calculation. In practice, it is computationally prohibitive to calculate the self-energy on the same fine k LMT O -mesh required for the LMTO part. Thus, a rather sophisticated procedure that includes several transformations of the self-energy Σ k nn ′ between different basis sets has been developed to interpolate the self-energy calculated by the GW part of the program on a coarse k GW -mesh to finer k LMT O -mesh used by the LMTO part of the program. 26 However, matrix elements of the self-energy (1) between states with high energy (≥ 2 Ry) often cannot be interpolated with sufficient accuracy from the k GW -mesh to the k LMT O -mesh 26 . Note that a fine k LMT O -mesh is required when describing metals. The latter is in part due to the long range of the LMTO basis set (e.g. the smallest eigenvalue of the overlap matrix can be of the order of 10 −10 ). In order to overcome this k-interpolation problem, the high-energy part (ε LDA kñ , ε LDA km > E xccut ) of the difference between the self-energy and LDA exchangecorrelation potential, was substituted with a diagonal matrix with the diagonal elements given by linear function of the LDA energy, △V xc nñ = a + b × ε LDA kñ . Here "∼" over the subscript denotes that the function is represented in the basis of eigenfunctions ψ LDA kñ of the LDA Hamiltonian (the LDA basis) with eigenvalues ε LDA kñ . The energy cutoff parameter E xccut is typically of the order of 2-3 Ry. The constants a and b are fitted from calculated △V xc nñ at lower energies. The results for the calculated quasiparticle (QP) energies usually depend weakly on the cutoff parameter E xccut or constants a and b. More details on the k-interpolation procedure and the method used to control its accuracy for bulk calculations could be found in Ref. [26], section II, subsection G.
In the LDA basis the optimum Hamiltonian H 0 (k) (2) reads After modifications of the matrix △V xc nm as outlined above, the Hamiltonian H 0 nm (k) in a form of Eq. (4) is used in the LMTO part of the program to obtain the QP wave functions and energies.
In the present work, we found that the QSGW method requires additional modification when the system exhibits extended regions with small electronic density, for example when modeling a metal/vacuum interface. Explicitly, the matrix elements △V xc nm in the LDA Hamiltonian of Eq. (4) leads to slightly improper mixing between occupied LDA states with energies ε LDA km ≤ ε F that are spatially concentrated in the metal region and "vacuum" LDA states with energies ε LDA kñ ≥ ε F + Φ extending to the entire volume of the system (here ε F is the Fermi energy). As a result, after diagonalization of the Hamiltonian (4), the QP occupied states have small tails that decay unphysically slow as a function of distance from the metal surface. Several reasons may contribute to this slow, non-exponential decay of occupied QP wave functions into vacuum, for example: remaining errors in the k-interpolation of the △V xc nm , non-completeness of the LDA basis, numerical errors, etc.
Development of a general procedure to construct optimum QSGW Hamiltonian for systems with vacuum regions in a way that guarantees correct exponential decay of occupied QP wave functions in vacuum is beyond the scope of this paper. However, in application to the specific case of Al, we can overcome this problem by straightforward truncation of the unphysical non-diagonal matrix elements △V xc nm with ε LDA km ≤ ε F and ε LDA kñ ≥ ε F + Φ, using the fact that the wave functions of bulk Al are rather well described by LDA (see Fig. 1). Specifically, we truncate all non-diagonal matrix elements of △V xc nm if the energies ε LDA kñ and ε LDA km satisfy following conditions: Here E c is a cutoff parameter. The idea of the method is to allow the occupied LDA statesm to be mixed by matrix △V xc nm only with unoccupied LDA statesñ with energy less then ε F +E c . Such procedure will prevent occupied states from mixing with extended vacuum states if E c < Φ [at least in the first order of the perturbation theory, if consider △V xc nm in Eq. (4) as a perturbation to the LDA Hamiltonian]. Note that second condition in (5) always allows the occupied states to be mixed between themselves. In the limit E c → ∞ there is no modification of the matrix △V xc nm and the method reduces to the standard QSGW approach. In the opposite limit, E c → +0, the method becomes an "unoccupied states eigenvalue-only" self-consistent GW method. The "unoccupied states eigenvalue-only" means that unoccupied subblock of the matrix △V xc nm is diagonal, so the unoccupied QP wave functions are always equal to the LDA wave functions and only QP energies are modified due to the diagonal matrix elements △V xc nñ . Thus, parameter 0 < E c < ∞ smoothly interpolates between these two methods.

III. RESULTS AND DISCUSSION
We begin with analyzing different characteristics of bulk Al and Al surfaces as function of parameter E c . Figure 1 shows the density of states (DOS) of the bulk Al calculated by LDA, full QSGW (that corresponds to the limit E c → ∞), modified QSGW as specified by Eq. (5) with small cutoff parameter E c = 0.1 eV, and the standard "eigenvalue-only" self-consistent GW methods. (In the "eigenvalue-only" self-consistent GW method all non-diagonal elements of the GW addition to the LDA Hamiltonian are neglected △V xc nm = △V xc nñ δñm, so the QP wave functions are always equal to the LDA ones.) It is evident that all three modifications of the GW method produce very similar DOS, somewhat different from the LDA result. The "eigenvalue-only" DOS is very close to the full QSGW DOS with minor deviations in the energy range from ε F − 1 eV to ε F + 3 eV. More importantly, the DOS obtained in the full QSGW (solid line) and DOS obtained by modified QSGW with small cutoff E c = 0.1 eV (dotted line) are almost indistinguishable on the figure. This means that the truncation (5) of the matrix elements △V xc nm with arbitrary cutoff parameter E c ≥ 0.1 eV practically does not change the bulk Al electronic structure. This is an important result suggesting that we can safely neglect erroneous non-diagonal elements of △V xc nm with ε LDA km ≤ ε F and ε LDA kñ ≥ ε F + Φ in surface calculations.
Next, we turn to the Al/vacuum interface systems. Figure 2 shows the electron density averaged over the x and y directions for an Al(111) surface using five dif- ferent values of E c . The calculations were performed for 4 Al and 8 vacuum layers (N A = 4 and N V = 8). The expected well behaved exponential decrease of electron density away from the metal surface is seen for E c = 1.36 eV, 2.72 eV, and 4.08 eV. Importantly, the densities obtained using these three E c are indistinguishable. On the other hand, for E c above 4.08 eV, the density begins to deviate from the normal behavior; it sharply increases near the center of the vacuum. The density calculated with E c = 6.8 eV even increases, at some z, when the distance from the metal increases. Similar results, independence on E c and correct exponential behavior of density in vacuum for E c ≤ 4 eV, and unphysical behavior for E c ≥ 4 eV, is found for Al(100) and Al(110) surfaces. We note that the unphysical behavior of the electron density only occur for small absolute density values, 4-5 orders of magnitude smaller then the density in metal region. Also, the QP energy bands depend only weakly on E c : even the energy bands calculated with parameters E c = 2.72 eV and 6.80 eV (not shown) almost coincide with each other for states with energies less then ε F + 2 eV and begin to deviate slowly for higher energies. Increasing E c from 2.72 eV to 6.8 eV results in an up-shift of the QP bands with energy above ε F + 4 eV by a value ranged from 0 to 0.2 eV, depending on particular band. Figure 3 shows the calculated work function, Φ, for Al(111) as function of the cutoff parameter E c . The value of the work function does not depend on E c , up to E c ∼ 4 eV but changes above this threshold. Similar behavior -independence of the work function on E c for E c ≤ 4 eV and rapid change above E c ∼ 4 eV threshold, was obtained for Al(100) and Al(110) surfaces. The  Figures 1-3, we conclude that (1) The DOS of bulk Al does not depend on E c for E c ≥ 0.1 eV; (2) for all three Al surfaces the electron density is exponential in vacuum, does not depend on E c for E c ≤ 4 eV, and demonstrates unphysical behavior for E c ≥ 4 eV; and (3) for all three Al surfaces the value of the work function does not depend on E c for E c ≤ 4 eV, and sharply changes for E c ≥ 4 eV. Therefore, in the range 0.1 eV < E c < 4 eV, E c is large enough for △V xc nm to include all important matrix elements (at least at the level of bulk Al), and simultaneously small enough to not include △V xc nm that erroneously mix occupied LDA states with vacuum LDA states. Importantly, the work function in this range does not depend on E c and thus could be taken as the true QSGW value of the work function. Therefore, in all calculations presented below the E c parameter is fixed and set to E c = 2.72 eV.
The range of applicability of the method described by Eq. (5) is limited to materials (such as Al) for which LDA wave functions are adequate so the matrix elements of △V xc nm with ε LDA kñ − ε LDA km ≥ Φ can be neglected. For any given metal, this condition can be verified on bulk level without performing time consuming surface calculations. The method is not applicable to metals (such as d -electron Fe and Cu) for which the matrix elements △V xc nm with ε LDA kñ − ε LDA km > Φ play significant roles. As mentioned above, further efforts are required to develop an universal QSGW-derived method applicable to Fe, Cu and other metals for which simple truncation of the matrix elements (5) does not work. Figure 4 shows the variation of the calculated work function, Φ, as function of slab thickness N A . The calcu-lations were performed with LDA and QSGW for unrelaxed surfaces using the following parameters: E c = 2.72 eV, N V = 6, (22 × 22) k LMT O -mesh, and (6 × 6) k GWmesh in the surface BZ (for the more anisotropic Al(110) surface we used (22×16) k LMT O -mesh, and (6×4) k GWmesh). For the LDA calculations we used the Barth-Hedin 7 functional. Φ oscillates as the slab thickness increases. These oscillations are well known 32 and can be attributed to quantum-size effects (QSE). The positions of the local maximums of the LDA Φ at N A = 5, 8, and 11 for the Al(111) surface, and minimums at N A = 7 for the Al(100) surface and at N A = 5 and 12 for the Al(110) surface are in agreement with previous DFT calculations 32 . The work functions obtained by the QSGW method show similar QSE oscillations. We estimate the uncertainty in our calculated for N A = 12 [N A = 14 for Al(110)] values of the work function due to the QSE as ±0.03 eV, which is larger then uncertainties due to other computational parameters like the number of k points in the k GW -mesh.
Delerue et. al 34 and Freysoldt et. al 20 found sizable renormalization of the GW self energy in thin semiconductor films due to the image potential at the interface; this effect is as large as 0.2 eV for the band gap of Si slabs with a thickness below 3 nm 34 . This is a finite size effect, different from QSE. On the other hand, for metallic films the image potential inside the metal slab is well screened, so it has only a minor effect on occupied states concentrated within the slab. Since the work function is mostly affected by occupied states, we do not expect a significant image potential induced correction to the value of the QSGW work function. Furthermore, because there is no image potential in the LDA approach, this assumption is supported by the similar behavior of Φ for QSGW and LDA as a function of slab thickness, see Fig. 4.
Results from our work function calculations for three Al surfaces are shown in Table I in comparison with experimental data and results from other theoretical studies. For all three surfaces, our LDA/CA results are relatively close to those obtained by other groups. The values of work function of Al(111) and Al(100) surfaces calculated by different groups using the GGA/BPE method also are relatively close to each other. Thus, one can conclude that the results obtained using a specific DFT functional are converged (within 0.1 eV or better accuracy) for different code implementations. On the other hand, Table I    When compared with data from experimental photoelectric measurements carried out under ultrahigh vacuum 43 , the work functions obtained using the QSGW method for Al(111), Al(100), and Al(110) surfaces differ by 0.07 eV, 0.05 eV, and 0.09 eV, respectively. All three differences are less then 0.1 eV and of the order of the sum of the theoretical and experimental error bars; we there-  Note that all calculated work functions presented in Table I (except G 0 W 0 [jel]) follow the increasing trend Φ(111) < Φ(110) < Φ(100), in agreement with data. This behavior is considered an anomaly; most other fcc metals instead follow Smoluchowski's rule Φ(110) < Φ(100) < Φ(111). The anomaly is caused by an increased p-atomic-like character of DOS at the Fermi energy in aluminum for the three surfaces, a behavior different to that of most fcc metals 36 . We in this context note that earlier experiments 44 reported 0.21 eV and 0.22 eV smaller values for Al(100) and Al(110) work functions compare to that of Ref. [43]. However, Grepstad et al. 43 suggested that this discrepancy could be due to higher impurity concentration, in particular oxygen, in the earlier experiment. Table I shows that Al work functions calculated using the G 0 W 0 method differ little (0.01-0.02 eV) from the converged QSGW results. The G 0 W 0 results presented in Table I correspond to using LDA/BH as starting point for GW iterations. Similar 0.01-0.02 eV deviations from the converged QSGW results were found for G 0 W 0 when instead using LDA/CA. We also note that the convergence of the GW iterations is fast, meaning that the ini-tial LDA wave functions are close to the converged QP wave functions; a conclusion supported by the similarity of the QSGW and QSGW(e-only) DOS for bulk Al shown in Fig. 1. The substantial differences to previous G 0 W 0 calculations by Morris et al 24 (G 0 W 0 [jel] line in Table  I) Table I) should therefore neither be attributed to the non-selfconsistency of the G 0 W 0 method nor to errors associated with the choice of particular DFT functional. Instead, we propose that the differences are due to the jellium approximation employed in both those studies 24,25 .
For more correlated materials such as Fe or Cu (where LDA and GW wave functions overlap less than they do in Al) we expect larger deviations of the work functions calculated by G 0 W 0 and QSGW methods as well as a stronger dependence of the G 0 W 0 results on the DFT functional used to calculate G 0 and W 0 .

IV. SUMMARY
We have applied the QSGW and G 0 W 0 methods to calculate the work functions of Al(111), Al(100), and Al(110) surfaces. The G 0 W 0 results differ from converged QSGW results by less then 0.02 eV and this small difference can be attributed to significant overlap of the LDA and QP wave functions. The QSGW results are in excellent agreement with experimental data taken under ultrahigh vacuum conditions. The calculated values of the work functions do not depend on the DFT functional used for the initial Hamiltonian H 0 . These results suggest that QSGW method can be used for reliable and accurate calculation of the work functions with accuracy of the order of 0.1 eV or better.
We found that modifications of the original QSGW method 10,11,26 are required in order to apply the method to the metal/vacuum surface. In particular, special care should be taken to control the errors originated form (slight) improper mixing of the occupied and vacuum states. In some simple cases, such as Al, where LDA wave functions are already a good approximation to the QP wave functions, simple truncation of corresponding matrix elements [see Eq. (5)] are enough to control these errors.
The truncation method is not applicable to metals such as d -electron Fe and Cu for which the matrix elements △V xc nm with ε LDA kñ − ε LDA km > Φ play significant roles. For any given metal, this condition can be verified by studying the bulk electronic structure, thus without performing time consuming surface calculations. Further efforts are required to develop an universal QSGW-derived method applicable to Fe, Cu and other metals for which simple truncation of the matrix elements is not adequate.