{'fraction_non_alphanumeric': 0.027343301847106965, 'fraction_numerical': 0.01870052392060882, 'mean_word_length': 4.573316283034953, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimized Effective Potential Using The Hylleraas Variational Method

In electronic structure calculations the optimized effective potential (OEP) is a method that treats exchange interactions exactly using a local potential within density-functional theory (DFT). We present a method using density functional perturbation theory combined with the Hylleraas variational method for finding the OEP by direct minimization which avoids any sum over unoccupied states. The method has been implemented within the plane-wave, pseudopotential formalism. Band structures for zinc blende semiconductors Si, Ge C, GaAs, CdTe and ZnSe, wurtzite semiconductors InN, GaN and ZnO and the rocksalt insulators CaO and NaCl have been calculated using the OEP and compared to calculations using the local density approximation (LDA), a selection of generalized gradient approximations (GGAs) and Hartree-Fock (HF) functionals. The band gaps found with the OEP improve on the calculated results for the LDA, GGAs or HF, with calculated values of 1.16eV for Si, 3.32eV for GaN and 3.48eV for ZnO. The OEP energies of semi-core d-states are also greatly improved compared to LDA.


I. INTRODUCTION
Density functional theory (DFT) is one of the great successes in tackling the many-electron problem 1-3 . Density-dependent exchange-correlation potentials, such as the local density approximation (LDA) or generalized gradient approximations (GGAs), are useful and surprisingly accurate for many properties of materials despite the known deficiencies of these functionals. Generally, the LDA slightly underestimates lattice constants, while GGAs slightly overestimate them and both underestimate single-particle band gaps in solids. [4][5][6][7] The exchange energy can be calculated exactly using the nonlocal Hartree-Fock (HF) formalism, rather than being approximated as part of the LDA or GGA. However pure HF grossly overestimates band gaps in solids and is considerably more expensive computationally because of the fully non-local nature of the exchange operator.
The central theorem of density functional theory 1 demonstrates that a local potential V (r) is sufficient to completely define a many-electron system in the ground state. A local potential that includes the effect of the exact exchange interaction was first proposed by Sharp and Horton 8 and solved nearly a quarter of century later by Talman and Shadwick. 9 This potential is known as the optimized effective potential (OEP), because a local potential is created to optimally represent the non-local HF exchange potential. The total energy of the OEP system is variational with the potential in the same way the LDA or GGA total energy is variational in the electronic density. This allows the construction of equations to describe the potential within a Kohn-Sham DFT formalism. The resultant single particle excitation energies agree much better with experiment than those calculated using the LDA, GGAs and HF. 3 The improved OEP description of excited states is a consequence of the freedom from self-interaction error of both occupied and unoccupied Kohn-Sham states 10 . Within the LDA and GGAs the potential experienced by all states includes some degree of self-interaction 11 . In HF the exchange term cancels the self-interaction error in the occupied states, but no such cancellation is present for the unoccupied states. This results in too small an electronic band gap using the LDA and GGAs and too large in HF 12-14 . There is a distinction between the fundamental band gap of a material and the corresponding Kohn-Sham gap. The Kohn-Sham gap (E KS g ) is defined as the difference between the eigenvalues of the highest occupied orbital and the lowest unoccupied orbital where i is the eigenvalue of the ith orbital and N is the total number of electrons in the system. Whereas the fundamental band gap (E g ) is defined as the difference in the ionization potential and the electron affinity where I is the ionization potential, A is the electron affinity and E(N ) is the total energy of a N electron system. The two are related by where ∆ xc is the derivative discontinuity in XC energy with respect to particle number. 6,7,15 This constant is the energy associated with a change in the exchangecorrelation potential when a infinitesimal charge is added to a N electron system. While a comparison between the Kohn-Sham and fundamental band gaps would be instructive, the value of the derivative discontinuity is unknown for most real materials. (See however Godby et al 16 and Gorling et al. 17,18 ) In finite systems, the exchange potential should decay as -1/r in the long-range limit for all states. In the OEP this decay occurs for all states irrespective of occupancy. In HF the non-local potential for occupied states behaves correctly but for the empty states decays exponentially and in the LDA/GGA the potential for all states decays exponentially 19 . Furthermore, the orbitals in the OEP each correctly decay with an individual exponent, whereas in HF all occupied orbitals decay with the same exponent. 3,10,[12][13][14] The OEP is potentially as versatile as LDA/GGA and HF methods. It could be used instead of HF as a foundation for orbital-dependent potentials (so called hybrid functionals) 3 . Furthermore once the local potential has been obtained, the calculation of system properties is faster than HF based methods, as the computationally expensive application of the exchange operator is avoided.
While the OEP correctly describes many system properties (see Grabo et al 12 ), it is not as widely used as the LDA, GGA and HF-hybrid based methods. This can be attributed in part to the lack of comparably efficient and robust computational schemes to evaluate it. In most formulations of the OEP method a sum over all excited states is required and truncation to a finite sum yields a slowly-convergent series. Furthermore an accurate description of these high energy states requires large basis sets at high computational expense if a local basis description is used. Despite these deficiencies, some calculations have been done in solids within the linear muffin-tin orbital (LMTO) basis set [20][21][22] , the FPLAPW method 63 as well as in plane-wave pseudopotential implementations 17,18,[23][24][25] . The agreement of the calculated results in the previous work above with experiment can be very good (see Ref. [10]). However most applications have been to semiconductors. For wide-gap insulating systems the performance of the OEP can be poor, notably for noble-gas solids. 26 See Ref. [10] for full details and further discussion of the OEP method and its results.
The extremely demanding task of calculating the full OEP can be reduced by using the mean-field approximation of Krieger, Li and Iafrate 27 (well-known as the KLI approximation) and even in this less precise formalism impressively accurate results have been obtained [23][24][25]28 .
The principal aim of this work is to derive and demonstrate a variational method of calculating the full exchange-only OEP without the need for a sum over all unoccupied states of a system. To this end, our method for the calculation of the OEP is formulated using density functional perturbation theory (DFPT) and the Hylleraas variational principle. As in the case phonon or electric field perturbations, only occupied Kohn-Sham orbitals are explicitly included. This method is then applied to a range of semiconductors and insulators.

A. The Optimized Effective Potential
We first summarize the OEP before deriving a variational implementation within the density functional perturbation formalism. For the remainder of this article we consider only the exchange-only OEP and will simply use the term OEP from here on. A starting point for finding the local exchange only OEP is the non-local Hartree-Fock equation whereT (r) is the kinetic energy,V ext (r) is the external potential,V H (r) is the Hartree potential,V X (r) is the HF exchange potential, i indexes the electronic states, σ is the spin index and φ σ i (r) is the ith orbital with spin σ. The Kohn-Sham (KS) equation 1,2 that incorporates the OEP is given by whereV σ is the effective potential fulfilling the role of both the Hartree and exchange potentials. The usual derivation of the OEP uses a chain rule expansion of the derivative of the exchange-correlation energy with respect to the density 12 , however here a treatment inspired by perturbation theory is used.
For a potential which differs from the ground state (GS) potential by an amount δV σ (r) such that V σ (r) = V σ (r) + δV σ (r) the GS KS orbitals change by where from first order perturbation theory 29 the last term of the above equation becomes The change in the effective potential also gives a first order change in the total energy, thus The OEP energy E[V σ ] is given by the Hartree-Fock functional evaluated using the Kohn-Sham OEP orbitals. Using the chain rule and the perturbation of the orbitals, the derivative of the energy with respect to the potential is given by By applying the Hellmann-Feynman theorem (for brevity the explicit r dependence of the potentials has been omitted) and substituting equations (7), (10) into equation (9), we obtain By substitution of the Kohn-Sham equation (eq. 5) the above can be written as 31 which is known as the OEP equation and was first derived by Sharp and Horton 8 . Here it is important to note that the requirement to sum over all unoccupied states presents a challenge for convergence of the system properties. A very large number of unoccupied orbitals must be included to ensure adequate convergence and consequently the calculation is extremely demanding 32 . The truncation of the infinite sum has also been the subject of concerns about the analytic properties of the solutions 33 .

B. Applying The Hylleraas Variational Principle
We present the following formalism to evaluate the OEP which avoids an explicit sum over states by borrowing ideas from density functional perturbation theory 29,37 using the Hylleraas variational method. The effective potential is obtained by variational minimization of E[V σ (r)]. The minimization direction for the potential is defined by the functional derivative of the energy with respect to the effective potential The functional derivative of the energy can alterna-tively be written as wherẽ is the first order correction to an unperturbed orbital φ σ i (x). These first order corrections to the orbitals were named "orbitals shifts" by Kümmel and Perdew. 13,34 An alternative method of calculatingφ σ i (x) without explicitly including any unoccupied states is to use the the Hylleraas variational principle 35 . We define a second order variational functional where˜ σ i is the first order correction to the eigenvalues.
and using the Hylleraas variational principle it can be shown that the orbital shifts which minimize the secondorder functional also satisfy equation (15). By substitution of the projection operator equation (16) may be rewritten The exact first order correction to the orbitalsφ σ i that minimizes G σ i [φ σ i ] can also be found using the Sternheimer-like equation 36 Equation (20) can be solved using iterative methods 37 . Replacing the infinite sum over states in equation (15) with an iterative procedure and considering only the occupied Kohn-Sham subspace avoids the severe convergence difficulties of the former method. Instead of determining of how many additional states to include in equation (15), at a computational cost increasing with their number, the convergence problem is transformed to one of determining how many iterative cycles are required.
Consequently a reliable solution may be achieved at considerably lower expense than summing equation (15) directly. By combining equations (14) and (20), the OEP can be found variationally while the first order correction to the orbitals is obtained by solving the Sternheimer equation directly.
Equations (14) and (20) are similar to those given by Kümmel and Perdew 13,34 . Our method differs from theirs by exploiting the variational character of E[V σ ] using the explicit gradient of equation (14) to perform a direct minimization. Their method involves a self-consistency cycle with an update procedure for v XC using a KLI-like expression and involving a division by the density. Our direct minimization is conceptually simpler and, because of the variational principle is likely to be more numerically robust.

III. IMPLEMENTATION AND CONVERGENCE
Our procedure for finding the OEP is to solve a double nested loop of minimizations. A flow diagram illustrating the full procedure is shown in Fig. 1. The inner of the two loops represents the solution of the Sternheimer equation (equation (20)) to find the first order correction to the orbitals. The method is derived from the Baroni Green's function solver technique, employing a conjugate gradient minimization scheme to findφ σ i (r) in a fixed potential V σ (r). [37][38][39] Using the first order correction to the orbitals and the orbitals themselves, the derivative of the total energy with respect to the potential can be found from Equation (14). This gradient is used in the direction-set methods to variationally optimize the potential, and its magnitude (residual norm) provides an indication of the level of convergence of the effective potential. The variation of the potential is accomplished via either one of two methods; the first is a conjugate gradient scheme 41 , the second method is a modified steepest descent method, known as the Barzilai-Borwein (BB) method, which requires knowledge of the gradient at the current point and previous point only 40 .
After each step of varying the effective potential, Kohn-Sham orbitals in this new, fixed potential are found non self-consistently. When the gradient and difference in total OEP Kohn-Sham energy per step is smaller than a predetermined threshold. the calculation is considered to be converged.
Calculate conjugate gradient search direction (g) for orbit Variationally solve eq. 20 fororbits in direction g The procedure is initialized by first solving the KS equation using a local density-dependent exchangecorrelation potential, in this case the LDA. From this self-consistent calculation the initial orbitals and density are obtained. Using this density the effective potential is then constructed explicitly using the LDA functional for the exchange and correlation. Solving the KS equation to find the ground state orbitals is accomplished variationally using a conjugate gradient method. 41 The non-local HF exchange energy and the Hartree energy are also initially calculated using the LDA orbitals.
The OEP has been implemented within the pseudopotential plane-wave code, CASTEP 42,43 . The orbitals, density and potentials are represented on rectilinear grids in the usual manner of a plane-wave DFT implementation. 45 Kohn-Sham orbitals are described on reciprocal space grid points G within a sphere bounded by the cut-off wavevector, G max and the density and potentials are nonzero within a sphere of radius 2G max . Therefore the grids in both real and reciprocal space used to represent the density and potentials have twice the dimensions of the grid used for the orbitals. The Hylleraas minimization scheme is performed explicitly on these real space grids by direct variation, so that the effective basis used to represent V σ (r) is the set of grid points {G} : |G| ≤ 2G max .
Optimized norm-conserving pseudopotentials used in this work were generated using the Opium code 46 developed by Rappe et al.. The Hartree-Fock approximation was used and the non-analytic behaviour of HF pseudopotentials was treated using the localization and optimization scheme of Al-Saidi, Walter and Rappe 47 . Including exchange exactly in the pseudopotentials is required to ensure accurate core-valence interaction within the current formalism. Incorrect core-valence interaction has been shown to have noticeable effects on calculated electronic structures. 48,49 We found that calculations using pseudopotentials which include d-states in the core typically predict band gaps 20-70% larger than if the dstates are treated as valence. Therefore semi core d-states were treated as valence for As, Cd, Ga, Ge, In, Se, Te and Zn.
The convergence characteristics of the conjugategradient and Barzilai-Borwein variational minimizers, are compared in Figure 2 which plots the total OEP energy for diamond as a function of iteration number. An initial rapid decrease in the total energy is followed by a long tail of decreasing change in the total energy in both cases. The conjugate gradient energy decreases smoothly and monotonically towards the converged result, but the Barzilai-Borwein energy does not, exhibiting sharp drops and an occasional spiked increase (behaviour which has been noted previously 50 ). After 100 or so steps the BB method still converges rapidly while the conjugate gradient method exhibits a very slow energy decrease. The difference in total energy per step is below 2.5µeV per atom after 250 iterations of the BB minimizer and the solution is stationary after 300 iterations. For the conju- gate gradients method however the energy difference per step does not drop below 5µeV per atom even after 500 iterations. After 500 iterations the energy difference between the two methods is within 8µeV per atom. In both cases we find that the calculated gradient tends to zero as the total energy converges, as expected in a variational method. The convergence rates of the Γ point band gap of diamond behave in a similar fashion as shown in Figure 3. After 250 iterations the gap is converged to within 2.5µeV per atom for the BB minimizer. By comparison, with the conjugate gradient minimizer the band gap is converged to 50µeV per atom in the same number of iterations. All calculations presented below using the OEP method were run for at least 250 iterations to ensure sufficient accuracy.
In the calculations that follow, the basis set size (planewave cut-off energy) and Brillouin-Zone sampling were chosen so that total energy differences, evaluated using the LDA, were less than 1 meV/atom. The same settings were used for the OEP calculations, which resulted in a very similar convergence error with cutoff as the LDA case. However the OEP calculations were more sensitive to Brillouin-Zone sampling error; for example, in the case of diamond, 8×8×8 resulted in a 1 meV/atom error using the LDA but 3 meV for OEP rising to 9 meV for ZnO. Convergence testing was also performed on the FFT grid used to represent the potential and the energy error was determined to be ≈ 10µeV.

IV. CALCULATED ELECTRONIC PROPERTIES
An example OEP is displayed in figure 4, in this case a slice through the primitive unit cell of Si in the diamond structure. The potential is smoothly varying outside the pseudopotential regions, which are the only regions where the potential is positive. The potential is at its deepest in the bonding region between the two ions and (outside the pseudopotential core region) is always negative but tending to zero in the low density regions away from the ions.
The electronic structure of a selection of insulating and semiconducting materials were calculated using the LDA, some GGAs (PBE 51 , PBESOL 51 , PW91 52 and WC 53 ), HF and the OEP. The resultant band gaps are shown in table I and for the LDA, OEP and Hartree-Fock calculations are also plotted in figure 5. The mean absolute error is 0.55eV for the OEP, 1.39eV for the LDA and 1.19eV for PBE, although these values are skewed by the underestimated values for the band gap found for diamond, NaCl and CaO materials. The corresponding band structures for Ge, CdTe, InN, GaN and ZnO are plotted in figures 6-10. Figure 6 shows the band structure of Ge in the diamond structure calculated using the LDA and the OEP. For Ge the OEP improves the gap from the LDA value of 0.09eV to 0.86eV which is very close to the experimental gap. Similarly for Si the gap is improved from 0.43eV using the LDA to 1.16eV using the OEP. For GaAs the OEP opens the gap from the LDA value of 0.99eV to 1.86eV. The OEP gap is greater than the experimental gap by 20%; a similar overestimation of the gap by the OEP in GaAs was also observed by Gorling et al. 17 Fig-ure 7 shows the band structures for CdTe calculated using the LDA and the OEP. The OEP gap is 35% larger than the experimental gap, opening from 1.46eV for the LDA to 2.20eV for the OEP. Figure 8 shows the band structures for the wurtzite semiconductor InN calculated using the LDA and the OEP. The calculated band gap for InN has a small direct gap of 0.21eV for the LDA which is opened to a gap of 1.39eV for the OEP. Similarly the gap in GaAs and CdTe is overestimated when compared to the experimental gap. The band structures calculated for GaN in the wurtzite structure are shown in figure 9. The OEP band gap of 3.32eV is very close to the experimental gap of 3.39eV and greatly improved on the LDA band gap 2.13eV. The wurtzite structure ZnO band structures calculated using the LDA and the OEP are shown in figure 10. The calculated band gap is greatly improved from 1.93eV when using the LDA to 3.48eV for the OEP which is very close to the experimental gap of 3.43eV.
For ZnSe in the zinc blende structure the OEP calculated gap of 2.83eV is remarkably close to the experimental gap and greatly improved over the LDA value of 1.88eV. However for diamond the OEP gap of 4.87eV underestimates the experimental value by 12% but still improves on the LDA value of 3.98eV, the underestimation by the OEP was again noted by Gorling et al. 17 For the wide gap rocksalt structure insulators CaO and NaCl the OEP band gaps are opened up compared to the LDA band gaps, with the LDA calculated band gap being 3.93eV for CaO and 4.84eV for NaCl and the OEP calculated band gap being 6.09eV for CaO and 6.32eV for NaCl. For CaO the gap is underestimated by 14% for the OEP and for NaCl the predicted gap is 30% smaller for the OEP than the experimental value. The OEP value for CaO is very similar to that found by Kotani and Akai 21 and for NaCl is in line with that found by Li et al. 64 Table II compares the OEP band gaps obtained here to those obtained by others using a variety of basis sets. The band gaps calculated here are very similar to previous values and confirm the progressive underestimation of the band gap for wide gap materials.
The OEP also improves values for semi-core energy levels relative to the LDA, as shown in table III. This is to be expected given that the self-interaction error present in the LDA has the largest effect on the tightly-bound d states, and that there is no self-interaction error within the OEP. In Ge the 3d electrons occupy states in the range -31.9eV to -32.0eV below the valence band maximum with the experimental binding energy being determined as 29.1eV to 29.6eV. The LDA puts them in the range -35.0eV to -35.2eV. The experimental In 4d electron binding energy in InN is 17.4eV; the LDA predicts between -18.4eV and -18.9eV and the OEP between -16.3eV and -17.2eV.
For GaAs the experimental Ga 3d states lie between 18.6eV and 19.0eV below the valence band maximum; the OEP predicts a range of between -20.0eV and -20.2eV and the LDA between -22.9eV and -23.0eV. Experiment      gives As 3d states at 40.4eV to 41.0eV below the valence band maximum; the OEP puts them at -44.5eV and the LDA predicts -47.5eV. In CdTe the Cd 4d electrons lie between, -10.2eV and -11.0eV for the OEP, -11.9eV and -12.3eV for the LDA, compared with an experimental value of 10.10eV below the valence band maximum. The experimental results also place the Te 5s electrons very close to the Cd 4d states with a binding energy of 9.2eV below the valence band maximum which the OEP also appears to predict. For ZnSe the Zn 3d electrons lie in the range -9.4eV to -9.7eV below the valence band maximum for the OEP, -11.6eV to -11.8eV for the LDA and 9.2eV from experiment. For the Se 3d electrons the LDA gives energies of -61.02 compared with -58.5 eV with the OEP and 55.5eV from experiment.
For GaN the LDA gives the Ga d-electrons lying between -21.2eV and -21.6eV, the OEP between -18.3eV and -18.9eV with an experimental value of 17.9eV below the valence band maximum. Using the OEP for ZnO gives Zn 3d states in the range -6.6eV to -7.6eV relative to the valence band maximum compared with an exper-imental value of 7.4eV and a range of -8.9eV to-9.7eV with the LDA. As can be seen from the band structures for the above materials in the valence, the s and p states are almost identical when using the LDA and OEP. It is the d states which display the greatest difference.
The OEP also improves upon the predicted electronic structure given by the GGA methods of PBE, PBESOL, PW91 and WC which underestimate the gaps of the materials investigated here and Hartree-Fock which greatly overestimates the gaps. Further work on magnetic metaloxides will be reported in a future publication.

V. CONCLUSIONS
A method of obtaining the OEP which treats the local exchange potential exactly without using a sum over all unoccupied states has been derived using the Hylleraas variational method and ideas borrowed from density functional perturbation theory. This allows for the calculation of the OEP using a variational minimization scheme in real space. The electronic structure of well known materials with a wide selection of band gaps have been calculated and the band gaps for semiconductors are found to be in good agreement with experimental values, although for the larger band gap materials, diamond, CaO and NaCl, the calculated band gaps are still underestimated by 10-30%. Hartree-Fock pseudopotentials were found to give more accurate results than LDA pseudopotentials. The absence of self-interaction error within the OEP is manifest in a better description of semi-core d-states compared to the LDA. Their energies with respect to the top valence band are much closer to experimental spectroscopic measurements than within the LDA.