Fast, efficient, and accurate dielectric screening using a local, real-space approach

We begin by reviewing the definitions of the polarizability and the dielectric response that screens an applied, external potential. We use Hartree atomic units throughout, such that the electron charge, electron mass, Planck’s constant and Coulomb’s constant are given by $e=m_{e}=\hbar=4\pi\epsilon_{0}=1$. We make use of the one-electron Green’s function

$\displaystyle g^{-1}(1,2)$

$\displaystyle=g_{0}^{-1}(1,2)-V(1,2)-\Sigma(1,2)$

(1)

where each numerical index denotes position, time, and spin, $g_{0}$ is the non-interacting Green’s function, $V$ is the total potential, which is local and so contains a factor $\delta(1-2)$, and $\Sigma$ is the self-energy encompassing many-body interactions.

The irreducible polarizability $\chi_{0}$ of the electron system is the change in electron density $n$ in response to a change in the total potential $V$,

$$\chi_{0}(1,2)=\delta n(1)/\delta V(2)\;.$$

(2)

The density can be written in terms of the one-electron Green’s function

$$n(1)=-ig(1,1^{+}),$$

(3)

where $1^{+}$ refers to a time infinitesimally later than $1$. Taking functional derivatives with respect to the potential in Eq. 1, the polarizability can be written

	$\displaystyle\chi_{0}(1,2)$	$\displaystyle=-i\delta g(1,1^{+})/\delta V(2)$
		$\displaystyle=\int\!\!d3d4\;ig(1,3)\frac{\delta g^{-1}(3,4)}{\delta V(2)}g(4,1^{+}).$		(4)

By approximating the Green’s function using $g^{-1}=g_{0}^{-1}-V$, we arrive at the random phase approximation for the polarizability,

$\displaystyle\chi_{0}^{\textrm{RPA}}(1,2)=-ig_{0}(1,2)g_{0}(2,1^{+}).$

(5)

This can be transformed from the two-time representation to the response as a function of a single external energy $\omega$ and written in real-space: $\chi_{0}(\mathbf{r},\mathbf{r}^{\prime},\omega)$.

Above, in Eq. 1, the potential term includes both the external and Hartree terms, and the Hartree term itself will change with changes in the electron density. One should therefore use the reducible polarizability $\chi$ which is the response to only changes in the external potential,

	$\displaystyle\chi(\mathbf{r},\mathbf{r}^{\prime},\omega)$	$\displaystyle=\chi_{0}(\mathbf{r},\mathbf{r^{\prime}},\omega)$
		$\displaystyle+\chi_{0}(\mathbf{r},\mathbf{x},\omega)v(\mathbf{x},\mathbf{x^{\prime}})\chi(\mathbf{x^{\prime}},\mathbf{r^{\prime}},\omega)$		(6)

where repeated spatial indices $\mathbf{x}$ are integrated over. Here $v$ is the Coulomb operator. Most importantly, from the reducible polarizability one obtains the screened Coulomb operator $W(\mathbf{r},\mathbf{r}^{\prime},\omega)$

	$\displaystyle W(\mathbf{r},\mathbf{r}^{\prime},\omega)$	$\displaystyle=\epsilon^{-1}(\mathbf{r},\mathbf{x},\omega)v_{\textit{ext}}(\mathbf{x},\mathbf{r}^{\prime})$		(7)
		$\displaystyle=v_{\textit{ext}}(\mathbf{r},\mathbf{r}^{\prime})+v(\mathbf{r},\mathbf{x})\chi(\mathbf{x},\mathbf{x}^{\prime},\omega)v_{\textit{ext}}(\mathbf{x}^{\prime},\mathbf{r}^{\prime}),$

where $\epsilon$ is the dielectric tensor. This is central to many-body perturbation techniques, for instance, for treating single-particle self-energies via the GW method or electron-hole excitation calculations via the BSE [12].

One can divide a potential into pieces and separately consider the screening of each piece. The two-coordinate external Coulomb operator in Eq. 7 can, without loss of generality, be written with a strength $q$ and parametric dependance on the second spatial coordinate

$$v_{\textit{ext}}^{[\mathbf{r}^{\prime}]}(\mathbf{x})=q|\mathbf{x}|^{-1}=q|\mathbf{r}-\mathbf{r}^{\prime}|^{-1}$$

(8)

where $\mathbf{x=r-r^{\prime}}$. Following Ref. [1], this potential is partitioned by adding and subtracting a shell of charge:

	$\displaystyle v_{\textit{ext}}(\mathbf{x})$	$\displaystyle=v_{1}(\mathbf{x})+v_{2}(\mathbf{x}),$
	$\displaystyle v_{1}(\mathbf{x})$	$\displaystyle=v_{\textit{ext}}(\mathbf{x})-v_{2}(\mathbf{x}),$		(9)
	$\displaystyle v_{2}(\mathbf{x})$	$\displaystyle=qR_{S}^{-1}\Theta[R_{S}-x]+qx^{-1}\Theta[x-R_{S}].$

Here $\Theta$ is the Heaviside theta function and $R_{S}$ is the shell radius. The screened potential is therefore

	$\displaystyle W^{[\mathbf{r}^{\prime}]}(\mathbf{r},\omega)$	$\displaystyle=\epsilon^{-1}(\mathbf{r,x=r-r^{\prime}},\omega)\,v^{[\mathbf{r}^{\prime}]}_{\textit{ext}}(\mathbf{x})$
		$\displaystyle=\epsilon^{-1}v_{1}+\epsilon^{-1}v_{2}$		(10)
		$\displaystyle=W^{(1)}\;+W^{(2)}.$

To this point no approximation has been made in the screening of $v_{\textrm{ext}}$. The full, two-coordinate $W(\mathbf{r,r^{\prime}})$ is recovered by evaluating the decomposed, real-space screening at each $\mathbf{r^{\prime}}$, as is shown for valence BSE calculations in section V.

The utility of the local, real-space approach rests on treating the screening of $v_{1}$ and $v_{2}$ differently. Because $v_{1}$ is nonzero only within $R_{S}$, the dielectric response that screens $v_{1}$ is also localized. By only calculating the polarizability within a finite volume, we are able to reduce the computational cost of the calculation. In contrast, the screening of $v_{2}$ must still cover all space. We therefore treat the screening of $v_{2}$ only approximately, using a model for the dielectric response.

Our method is therefore an approximation to the RPA screening, controlled by the shell radius $R_{S}$ and the model screening used to screen $v_{2}$. As $R_{S}$ increases, $v_{2}$ goes to zero, and treating $v_{2}$ approximately is a controlled approximation. The effect of dividing the calculation at finite $R_{S}$ is addressed further in section IV.3 and appendix D.3. A modified Levine-Louie dielectric model was adopted in Ref. [1] and is used here [13, 14]. This model screening is parametrized by the static, long-range dielectric constant $\epsilon_{\infty}$ and the average and local valence electron densities.

This choice makes $\epsilon_{\infty}$ an input parameter, but small errors in $1/\epsilon_{\infty}$ will only lead to small errors in the resulting x-ray absorption spectra, and we investigate this further in section IV.3 and appendix D.4. As suggested in Ref. [1], the preferred source is an experimental measurement of the optical constants or index of refraction below the band gap. If such data are not available, several computational approaches can be used. The optical spectra can be calculated with the ocean code [15, 16], either within the RPA or within the Bethe-Salpeter equation (BSE) approximation, which formally requires iteration to equate the input and output $\epsilon_{\infty}$.

The short-ranged potential $v_{1}$ is screened by calculating the RPA polarizability. Typically we use a shell radius $R_{S}$ between 3 a.u. and 5 a.u. The polarizability is calculated within a spherical region of space given by a radius of 8 a.u. to 10 a.u. Per Eq. 5, this necessitates calculating the one-electron Green’s functions for this region of space. Improvements in the fidelity, efficiency, and parallelization of calculating $g$, $\chi_{0}$, and $\chi$ in order to screen $v_{1}$ are the focus of this work.

The screening method of Ref. [1], adopted in the ocean spectroscopy code, utilizes electron wave functions generated from a pseudopotential-based density-functional theory code. The pseudopotential approximation allows for a dramatic reduction in the computational cost of planewave codes. The core-level electrons are removed, and the $-Z/r$ potential of each ion is replaced by a pseudopotential, such that the valence electrons are the most-bound states. This reduces the computational by treating fewer electrons and smoothing the remaining electron wave functions, reducing the required number of planewaves. More information on pseudopotential theory and methods can be found in [17].

Consider a pseudopotential $V^{\textit{ps}}_{\textrm{ion}}$ for an element which treat angular-momentum-dependent effects separably following the Kleinman-Bylander form [18]. For each value of the principle angular momentum $l$ up to $l_{\textrm{max}}$ it consists of some local potential and optionally some non-local projectors. In contrast, the all-electron ionic potential is simply an attractive Coulomb potential set by the ion’s atomic number $Z$. Both the all-electron ($ae$) and pseudized ($ps$) radial Schrödinger equations can be solved numerically.

		$\displaystyle H^{ae/ps}\phi_{j}^{ae/ps}=\varepsilon_{j}\phi_{j}^{ae/ps}$		(11)
		$\displaystyle H^{ae/ps}=-\frac{\nabla^{2}}{2}+\frac{l(l+1)}{2r^{2}}+V_{H}[n]+V_{xc}[n]+V^{ae/ps}_{\textrm{ion}}\,$

The Hartree $V_{H}$ and DFT exchange-correlation $V_{xc}$ potentials both depend on the density $n$, and so the problem must be solved self-consistently. Outside of some cutoff radius $r_{c}$ the pseudopotential matches the all-electron ionic potential, and the pseudo and all-electron electron orbitals are equal for an isolated atom. Additional requirements are enforced at $r_{c}$. The all-electron and pseudo orbitals should have the same radial derivative and scattering length, although this is exact only for specific energies. Because the scattering properties are the same for the all-electron and pseudo wave functions, the behavior of the pseudo electrons in the interstitial regions is identical to that of the electrons in the all-electron system.

While pseudopotentials are capable of reproducing all-electron results for many properties, such as band structures (Bloch-state energies) or structural properties, they fall short in others, such as core-level excitations. This is unsurprising because there are no core-level electrons in the pseudo-potential system. We point to localization and location of a perturbation or transition operator in determining if the pseudopotential wave functions are sufficient. In the case of core-level transitions, the operator is both highly-localized and within the pseudized region. Conversely, structural properties are not localized. They rely on force constants, which are a measure of the change in the distribution electron density as felt by one atom in response to the motion of another. In the case of a real-space screening approach, we are interested in the density response to a localized perturbation, and for screening a core hole this localized perturbation is within the pseudized region.

To demonstrate the differences between an all-electron and pseudopotential screened core-hole, in Fig. 1(a) we show the difference in the valence (2s and 2p) electron densities between the all-electron and the pseudopotential in a fluorine atom. While the all-electron 2s orbital has a node (and is orthogonal to the 1s orbital), the pseudopotential 2s does not. The first anti-node of the all-electron 2s is responsible for the difference in densities below around 0.2 a.u. Around 0.7 a.u. the pseudopotential density exceeds the all-electron density, a consequence of norm conservation and the deficiency at small radii. In Fig. 1(b) we plot the change in valence electron density upon introduction of a 1s core hole in the system. For both systems the valence electron density moves to smaller radii (towards the positively charged core-hole), but for the all-electron system there is a larger change at short distances.

The small difference in density response is magnified when converted into an induced potential as shown in Fig. 1(c). Clearly, at short distances the all-electron valence orbitals are more efficient at screening core-level excitations than the pseudopotential orbitals. In section IV we will explore how this difference in induced potential affects calculated x-ray absorption spectra that can be compared to measured ones.

In the context of screening a core hole, the potential $v_{\textit{ext}}$ is already weakened by relaxation of the core orbitals that are not included explicitly when valence screening effects are computed. This largely prevents problems that could arise from a frozen-core approximation regarding changes in the core level occupancies. Smaller, additional core-relaxation effects arise when there are changes in the chemical environment. These can include core polarization (mostly dipolar) and relaxational self-energies of valence electrons or holes added onto an atomic site. One method to treat many of these effects is presented by Shirley and Martin [19], who also refer to previous work.

To include the correct all-electron behavior, we augment wave functions of the pseudopotential system. This augmentation style was introduced for core-level transitions [20, 15] and borrows heavily from the projector-augmented wave (PAW) method [21]. In the PAW method, Bloch states are augmented by projecting smoothened wave functions onto a basis set of projector functions. For DFT ground-state calculations, it may be advantageous to optimize the projection scheme for describing the highest occupied and lowest unoccupied states, whereas our emphasis is obtaining realistic electron wave functions over a wide range of occupied and unoccupied bands. The coefficient of each projector function weighs the correction of the wave function in the form of replacing the projector function with an all-electron counterpart.

In ocean, all-electron and pseudo versions of partial waves are evaluated and condensed into a set of optical projector functions (OPFs). Typically, partial waves are sampled at several dozen regularly spaced energies over a multi-hartree energy range. The partial waves are constructed within an augmentation radius the encompasses the the pseudotpotential cut-off radius $r_{a}\geq r_{c}$. The OPFs are the eigenvectors that have the largest eigenvalues of the overlap matrix between all of the partial waves, considering each angular momentum separately. As a result, the OPFs are not single-energy partial waves, in contrast to typical PAW construction. The eigenvalues fall of rapidly enough that over 99.9 % of the wave function’s degrees of freedom are sampled with only a few OPFs. This is described in further detail in appendix C.

Augmentation exploits the same properties that we earlier asserted our pseudopotential would have. Namely, outside of the pseudopotential cut-off radius $r_{c}$ the pseudo wave function is identical to the all-electron wave function. Over a reasonable energy range, a bound or scattering state can be transformed from the pseudopotential system to the all-electron system by replacing the pseudo wave function inside of the cut-off radius with the all-electron solution. The OPF method is only used for augmentation to reconstruct wave functions as a means of post-processing DFT results, and the implementation is discussed in section VI.3.

Previously, the ocean code relied on an approximation instead of carrying out augmentation of the wave functions. We document the old approach here, and in the next section we will compare it to the current method. As was shown in Fig. 1(c), even the isolated atom demonstrates the importance of all-electron wave functions when calculating the screening near the nucleus. The difference between the two induced potentials in Fig. 1(c) can be calculated purely within the isolated atomic case

$$\Delta v_{\textrm{ind}}(r)=v^{ae}_{\textrm{ind}}(r)-v^{ps}_{\textrm{ind}}(r)\;.$$

(12)

Previously this correction was at times applied to the induced potential as calculated within the RPA using the un-augmented wave functions of the system. In this way the approximate effect of augmentation was included in the screening. However, because this method relies on calculations of an isolated atom, the valence orbitals and their occupations are only a rough approximation to the system of interest.

Before showing the effects of augmentation on calculated spectra, we examine the effects on the screening, using LiF as a model system. In Fig. 2 we show the spherically averaged induced potential in response to a F 1s core hole. As in the case of the isolated atom (Fig. 1), the induced potential is stronger near the origin when the orbitals have all-electron character, and the effect of augmentation is only significant near the atom where the all-electron and pseudo-potential orbitals are different. The approximate augmentation, while exact for an isolated atom, does not fully reproduce the screening using augmented orbitals within a solid. It is also only valid for atom-centered response such as for core-level spectroscopy.

For valence-level spectroscopy or other calculations, such as self-energy calculations using the GW method, the screening must be calculated throughout the unit cell. In Fig. 3 we show the induced potential in response to test charges centered along the $(a/2,a/2,a/2)$ (from F to Li) and $(0,a/4,a/4)$ (from F halfway to the nearest F) directions within the unit cell, where $a=4.0$ Å. Unsurprisingly, the effect of augmentation is only noticeable near the F atom. Note that while the point $(a/2,a/2,a/2)$ is centered on a Li atom, there is no effect from augmentation because the Li 1s orbitals are included as valence.

We first consider the fluorine K edge in LiF, which has been studied with ocean previously [15, 30]. Fig. 4 shows the effect of changes to the screened Coulomb attraction between the electron and core hole. In light grey, the non-interacting spectrum (neglecting electron-hole interactions) is dramatically different from any other approximation, and it shows the importance of the excitonic binding on the absorption spectrum. The three different approximations to the screening all give qualitatively the same results, capturing a strong exciton around 1.5 eV below the onset of the non-interacting spectrum (0 eV in the plot). The differences are mostly confined to the near-edge region, within 10 eV of the onset. Without any augmentation the exciton is 1.52 eV below the conduction band minimum. This reduces to 1.34 eV with approximate augmentation and to 1.29 eV with true augmentation. Correspondingly, the exciton strength is reduced with augmentation, improving agreement with experiment (Fig. 5).

Small changes in intensity and spacing of peaks in the x-ray absorption spectra can be signs of changes in the local structure around the absorbing atom [31]. To illustrate the importance of accurately determining the relative intensities of peaks, we reproduce measured XAS from reference [30] which shows the change in the F K-edge absorption with temperature. LiF has a dipole-forbidden pre-edge that is observable in dipole-limited absorption spectra due to the vibration of the ions which breaks the symmetry around the F atoms [30, 31]. When simulating spectroscopy of condensed systems, disorder is typically treated within the Born-Oppenheimer approximation where the ions are stationary during the excitation process. The final calculated spectrum is an average over spectra from different atomic positions, e.g., generated using molecular dynamics simulations [31] or statistical sampling of the vibrational modes of the system [32], or a statistical average of core-hole displacements [33]. Typically in calculations of extended system, the vibrational energy levels of the ionic system are treated as negligible and Frank-Condon effects are neglected [33].

In Fig. 5 we compare the calculated changes in the F K near-edge spectra of LiF due to changes in the screening calculation to measured changes due to temperature. At increased temperature, the pre-edge feature near 692 eV moves to lower energy, all of the features are broadened, and there is a noticeable weakening of the main exciton between 694 eV and 696 eV. The main exciton also differs in strength between the calculations using different augmentation methods, which is even more exaggerated if augmentation is neglected (Fig 4). While the error in the calculation from neglecting or approximating the augmentation has only a small qualitative effect on the spectrum, the differences are substantial when compared to small structural changes. High-fidelity screening calculations are necessary for correctly identifying local structure or assessing the accuracy of approximations used for incorporating disorder or vibrations.

Having shown the large effect of changes to the screening on the XAS of LiF, we next examine the effect of augmentation on heavier ions by exploring the range of lithium halides. All crystallize in the rocksalt Fm$\bar{3}$m structure, and for all four materials a uniform core-hole lifetime broadening of 0.5 eV was used. In Fig. 6 we show the effects of different augmentation approaches on the halide K edges of LiCl, LiBr, and LiI. The trend between approximations for the same compound shown for LiF in Fig. 4 holds for the heavier halides as well. Calculating the screening of the core-hole potential with wave functions from a pseudopotential calculation dramatically under-screens the core hole, resulting in an exaggerated excitonic peak. There is a trend towards smaller discrepancy with increasing atomic number in the halide series. For the LiCl the exciton is approximately 0.56 eV (0.26 eV) more bound without any augmentation (with approximate augmentation), while for LiBr the over-binding is 0.65 eV (0.11 eV) and 0.05 eV (0.01 eV) for LiI. We conclude that some all-electron augmentation is necessary for the proper calculation of the screening even for heavier atoms, but for Br and I, it appears that the approximate augmentation method may be sufficient.

The differences are primarily confined to the near-edge region, which could be expected from Fig. 1(c). As for LiF, the differences between the augmented and un-augmented induced potential are confined to a small region around the core hole. Only near the edge onset is the excited electron localized enough to be strongly affected by this very localized difference in core-hole potential. However, an increase in spectral weight near threshold necessitates a reduction in spectral weight higher in energy – high energy states in the non-interacting system are pulled down by the core hole. The differences shown here would be mitigated somewhat by more realistic core-hole lifetime broadening for the Br and I K edges, though by observing partial fluorescence emission from the $2{\mathrm{p}}_{\nicefrac{{3}}{{2}}}$, the 0.5 eV broadening used here remains realistic for the Br.

Next, we examine the effects of augmentation by comparing calculations of LiI using two different iodine pseudopotentials. The standard iodine pseudopotentials uses a Kr core with the 4d, 5s, and 5p electrons in the valence bands while the semi-core pseudopotential also includes the 4s and 4p orbitals as valence. In Fig. 7 we show that the calculated I K edge of LiI does not depend on the pseudopotential when the orbitals are properly augmented for the screening calculation. However, without augmentation the screening calculated using the standard pseudopotential is notably weaker than that of the semi-core pseudopotential, leading to a stronger exciton. The weaker effect of augmentation on the calculation using the semi-core pseudopotential is due to the node in the 5s and 5p orbitals (absent for the standard pseudopotential). As discussed earlier (see Fig. 1), the absence of nodes in the pseudopotential orbitals shifts the electron density away from the atom. Even though without augmentation the $n=5$ orbitals only get a single node in the semi-core system the effect is already dramatic and the discrepancy due to neglecting augmentation is strongly reduced.

We have demonstrated that the pseudopotential-based orbitals are insufficiently accurate for calculating the screened electron-hole interaction for core-level spectra. This failure is independent of the quality of the pseudopotential used. Instead it is a straightforward consequence of the difference in nodal structure between valence orbitals in the all-electron and pseudopotential systems. As shown in the atomic case (Fig. 1), the absence of nodes in valence orbitals of the pseudopotential system affects the density response and hence the screening. This is mitigated somewhat by the use of semi-core pseudopotentials that include the next highest principle quantum number in the valence, e.g., Fig. 7. The moderate success of including semi-core orbitals without augmentation may indicate that a kind of false augmentation could be used, where nodes are added without requiring that they mimic the true all-electron orbitals. These false projectors could be constructed such that they add nodes without increasing the required plane wave energy cutoff for use in reciprocal-space methods.

We next investigate the K edges of hexagonal boron nitride (h-BN). The XAS of h-BN has been the subject of recent investigations into the role that vibrational disorder and defects play in the spectra [34, 35, 36]. Because h-BN is layered, it has an anisotropic dielectric response, unlike the lithium halides. We will use it to showcase the real-space screening method on systems with reduced symmetry. The dielectric screening within the BN planes (perpendicular to the c-axis) is stronger than that parallel to the c-axis: $\epsilon_{\infty}^{\perp}=4.95$ versus $\epsilon_{\infty}^{\parallel}=4.10$ [37]. The anisotropy is also visible in the XAS. Both the B and N edges show strong excitonic features when the x-ray polarization vector is aligned with the c-axis and a delayed onset when the x-ray polarization is in-plane.

We compare our calculations of h-BN to two different computational approaches. First, the exciting code (like ocean) uses the BSE [38, 39]. The screening in exciting is also carried out in the RPA, but calculated in reciprocal space with no modeled response. Additionally, it is an all-electron code making augmentation unnecessary. Second, within OptaDOS x-ray spectra are calculated using the $\Delta$SCF method [40]. In this approach the final states are calculated directly as the unoccupied states of a DFT calculation with a core hole. The initial state is a the core-level orbital, and the transition matrix elements can be calculated directly. The density response to the core hole is calculated self-consistently within DFT. The OptaDOS calculations reproduced here used the castep DFT code [41].

In Fig. 8 we compare the N K-edge of h-BN calculated with ocean to experiment [42] and calculations using OptaDOS and exciting. The experiment was taken with polarization at an angle of $50^{\circ}$ to the c-axis, giving nearly the same 2:1 in-plane to out-of-plane ratio of a disordered sample. All three calculations are averaged over polarization directions. The $\Delta$SCF OptaDOS spectrum is a clear outlier with a substantially reduced exciton strength and binding at the onset, a known characteristic of $\Delta$SCF methods [44]. The agreement between ocean and exciting is very good. The exciting calculation includes fewer conduction bands, leading to an artificial die-off in spectral intensity at higher energies.

Per section II.3, the real-space screening method is only an approximation to the RPA response. The approximation is controlled by $R_{S}$, the radius of the shell of charge that is screened with a model dielectric response, and $\epsilon_{\infty}$, the value of the static dielectric constant which is input into the Levine-Louie dielectric model (see appendix B). In Fig. 9(a) we show the total induced potential in h-BN due to a core hole on the nitrogen atom. We also show the difference in the potential calculated using a shell radius of 7 a.u. versus smaller radii, and we find that the errors are less than 3 mHa. Next we vary the value of the input dielectric constant by approximately $\pm 15\%$, plotting the difference in the potentials in Fig. 9(b). As expected, increasing the dielectric constant input into the model decreases the induced potential. However, we find the differences are $\lesssim 5$ mHa. The screening calculation is relatively insensitive to changes in the sphere radius or errors in the input dielectric constant. Further examination is presented in appendix D.

The primary effect of small errors in the screening is a shift of the spectra, reflecting a change in the excitonic binding energy. For the nitrogen K edge, using $\epsilon_{\infty}=5.4$ blue-shifts the x-ray absorption by 0.10 eV, while using $\epsilon_{\infty}=4.0$ red-shifts the x-ray absorption by 0.12 eV. In Fig. 9(c) we show how the nitrogen K-edge spectra change with these changes in the calculation of the screening. The changes are minor, and so only the differences between the spectra are plotted, after taking into account the aforementioned red or blue shifts.

Previously in ocean, the screened potential for valence calculations was approximated using the Hybertsen-Levine-Louie dielectric model [13, 48], which depends parametrically on the local density $\rho(\mathbf{x})$ and static dielectric constant $\epsilon_{\infty}$. The wave functions for the electrons and holes are sampled on regular real-space grids $\mathbf{x}$, and therefore we need a description of the screened Coulomb potential for each set of grid points $W(\mathbf{x},\mathbf{x}^{\prime})$. Following Ref. [45], the screened Coulomb potential is given by,

	$\displaystyle W_{\textrm{HLL}}(\mathbf{x},\mathbf{x}^{\prime})=\nicefrac{{1}}{{2}}$	$\displaystyle\left[W_{\textrm{hom}}(|\mathbf{x}-\mathbf{x}^{\prime}|;\rho(\mathbf{x}),\epsilon_{\infty})\right.$
		$\displaystyle\left.W_{\textrm{hom}}(|\mathbf{x}-\mathbf{x}^{\prime}|;\rho(\mathbf{x}^{\prime}),\epsilon_{\infty})\right]$		(13)

which simply averages the results using the density at points $\mathbf{x}$ and $\mathbf{x}^{\prime}$. To avoid the divergence at $\mathbf{x}\rightarrow\mathbf{x}^{\prime}$, a spherical average over the discretization volume is used when $\mathbf{x}=\mathbf{x}^{\prime}$.

To improve this, we substitute the more accurate local-RPA result for the short-range part of $W$. Using the previously introduced method, we can calculate the screened Coulomb from Eq. II.2 for each grid point $\mathbf{x}$, ie, $W^{[\mathbf{x}]}$. In the case of core-level excitations the external potential in Eq. 7 was given by the core-hole potential from an atomic calculation. For valence calculations we use the potential from a spherical charge centered at $\mathbf{x}$. The volume of the spherical charge is set by the discretization volume of the real-space grid, $V_{x}=\Omega/N_{x}$, the unit cell volume $\Omega$ divided by the number of grid points $N_{x}$. The screening calculation can be carried out with or without augmentation.

As above, we enforce the symmetry in interchanging $\mathbf{x}$ and $\mathbf{x}^{\prime}$.

$$W(\mathbf{x},\mathbf{x}^{\prime})=\begin{cases}\frac{1}{2}\left[W^{[\mathbf{x}]}(\mathbf{x}^{\prime})+W^{[\mathbf{x^{\prime}}]}(\mathbf{x})\right],&\!\!\text{if }|\mathbf{x}-\mathbf{x}^{\prime}|\leq r_{m}\\ W_{\textrm{HLL}}(\mathbf{x},\mathbf{x}^{\prime}),&\text{otherwise}\end{cases}$$

(14)

As for the core-level, the RPA screening is calculated only within a finite space and parametrized by the shell radius $R_{S}$ and dielectric constant $\epsilon_{\infty}$. At large distances the HLL approximation is still used, governed by the parameter $r_{m}$. Much like the shell radius, this constitutes a controlled approximation.

Like the HLL model, care must be taken when evaluating the real-space $W$ in the limit of $\mathbf{x}\rightarrow\mathbf{x}^{\prime}$. We numerically integrate the $l=0$ component of the calculated screened potential over the discretization volume

$$W^{[\mathbf{x}]}(\mathbf{x})=\frac{3}{R_{x}^{3}}\int_{0}^{R_{x}}W^{[\mathbf{x}]}(r)r^{2}dr\;,$$

(15)

where $R_{x}=[3V_{x}/4\pi]^{1/3}$.

The system-size scaling behavior is the same for the valence screening as it was for the core-level case — the number of grid points $\mathbf{x}$ scales linearly with volume the same as the expected number of atomic sites. (The scaling is discussed in section VI and illustrated in VI.7.) There are two major differences between calculating the screening for the valence and core cases with negative and positive impacts on run time. First, the number of real-space grid points is much larger than the number of atoms. For example, in a unit cell of LiF an $10^{3}$ x-point mesh is necessary, resulting in 1000 screening sites instead of the 1 needed for the fluorine x-ray absorption calculation. The dramatic increase in the number of sites is offset by the use of coarser real-space grids in the calculation of $\chi$. The perturbing potential for the core case is the core-hole potential, and, like the core-hole density, it is strongly localized. In the valence case, the perturbing potential is taken to be a uniform ball of charge whose volume is set by the discretization volume $V_{x}$ defined above. Therefore, the valence screening calculations converge with a coarser radial mesh than is needed for the core-level screening.

To show the effects of using the local RPA, with and without augmentation, instead of the HLL screening on valence calculations we consider calculations of the imaginary part of the dielectric function.

First, we look at bulk silicon — a standard for valence electronic structure calculations providing a testbed for early DFT, GW, and valence BSE calculations [51]. As for the lithium halides we use the experimental lattice constant of 0.543 nm [52] and the PseudoDojo pseudopotential for silicon. Additional input parameters are included in appendix A. A scissor correction was used to set the DFT band gap to be 1.11 eV [53]. Unsurprisingly, we see in Fig. 10 that in the case of silicon the HLL model performs very well. The inclusion of the RPA screening at short range has little effect on the spectrum of bulk silicon, and augmentation has no visible effect.

Next we return to LiF (also well studied previously within the BSE [54, 15, 55]). The simulation details are similar to the x-ray case, but we have included a scissor correction to set the DFT band gap to 14.2 eV [56]. We first compare our calculations BSE results calculated using ABINIT v. 8.10.2 [57] and yambo v. 4.5.0 [58] codes in Fig. 11(a). A small, 0.09 eV shift in the position of the first exciton between the ocean and yambo indicates slightly weaker screening within ocean. The ABINIT results were obtained with a scissor-corrected band gap of only 13.3 eV, chosen to align the exciton position with the yambo calculation. The origin of the 1 eV shift in the ABINIT spectrum is not known. The three calculations are largely in agreement with each other in shape and intensity as well as previously published calculations, accounting for the choice of scissor correction.

In contrast to bulk silicon, LiF is strongly ionic with the valence orbitals primarily localized on the fluorine site, making it a more difficult system for the model dielectric response. In Fig. 11(b), we see discrepancies between the HLL-only and RPA results. Both the main exciton near 12 eV and the higher energy peak near 21 eV are red-shifted in the HLL as compared to the RPA calculations. Using the HLL model, the main exciton is at 11.8 eV, while the RPA places it at 12.0 eV. This shift indicates that the HLL model screening is weaker than the RPA calculation for these two peaks. This is consistent with a larger $GW$ band-gap correction that often occurs with the HLL approximation is used as compared to when the RPA is used [59]. As for silicon, the differences in the spectra with and without augmentation are very minor and not distinguishable in the plot.

The initial step is determining the electron wave functions and the basis, from which we can generate the Green’s functions. Density-functional theory (DFT) is used to simulate the electronic Hamiltonian. The system is taken to be periodic such that the electron orbitals can be denoted by their band $b$, crystal momentum $\mathbf{k}$, and spin $\sigma$,

$$H^{\textrm{DFT}}\psi_{b\mathbf{k}\sigma}=\psi_{b\mathbf{k}\sigma}\varepsilon_{b\mathbf{k}\sigma},$$

(16)

where each wave function $\psi_{b\mathbf{k}\sigma}$ has energy $\varepsilon_{b\mathbf{k}\sigma}$. The Green’s function for energy $E$ can be written in the spectral representation

$$g_{\sigma}(\mathbf{r},\mathbf{r}^{\prime},E)=\sum_{b=1}^{\infty}\int_{BZ}\frac{d^{3}\mathbf{k}}{(2\pi)^{3}}\,\frac{\psi^{\dagger}_{b\mathbf{k}\sigma}(\mathbf{r})\psi_{b\mathbf{k}\sigma}(\mathbf{r}^{\prime})}{E-\varepsilon_{b\mathbf{k}\sigma}}$$

(17)

The integral runs over the Brillouin zone.

To construct $g$ we must define the real-space basis. As mentioned in Sec. II.3, we calculated the RPA response only for the local potential $v_{1}$. We employ a real-space basis within a sphere $S$ with a radius $r_{S}$ and centered on a point $\tau$. For screening the core hole, $\tau$ is the atomic site, while for valence calculations $\tau$ is one of the grid points in Eq. 14. The irreducible polarizability $\chi_{0}^{(\tau)}$ is then an $N_{r}\times N_{r}$ size matrix, as are the Green’s functions $g^{(\tau)}(E)$. This real-space basis is independent of the size of the system’s unit cell, and is discussed further in Appendix D.1.1.

In practice the sum over bands is truncated, and the integral is replaced with a sum over regularly spaced points in k-space. Our approach requires only a few k-points, often between 1 and 8, and we address this later in Appendix D.2 on errors and convergence. The sum over bands, however, is significant. Typically, convergence in the screening is reached when the Green’s function is constructed with states up to around 100 eV above the Fermi level. This requires on the order of 30 to 50 bands per atom, and the number of bands required scales linearly with system size. Unfortunately, aspects of the generation of the one-electron states from DFT scale with the square of the number of bands.

These DFT eigenstates $\psi$ are generated using an external plane-wave DFT program and are saved to file as Bloch states $u$,

$$\psi_{b\mathbf{k}}(\mathbf{r})=e^{i\mathbf{k}\cdot\mathbf{r}}u_{b\mathbf{k}}(\mathbf{r})=e^{i\mathbf{k}\cdot\mathbf{r}}\sum_{\mathbf{G}}C_{b\mathbf{k}}(\mathbf{G})e^{i\mathbf{G}\cdot\mathbf{r}}$$

(18)

which are defined in terms of complex-valued coefficients $C$ of plane waves $\mathbf{G}$. The spin index $\sigma$ will be dropped. Only the set of coefficients $C$, not the various phases, are written to file. We distribute the work for the conversion by band and k-point — not by plane wave coefficient.

To project the wave functions onto our spherical grid, one option would be to follow the method given by Eq. 18 directly, i.e., Fourier interpolation. We first create a matrix of the phases, as these will be common across all of the bands at a specific k-point.

	$\displaystyle P_{\mathbf{k}}(\mathbf{r},\mathbf{G})$	$\displaystyle=e^{i(\mathbf{k+G})\cdot\mathbf{r}}$		(19)
	$\displaystyle\psi_{\mathbf{k}}(\mathbf{r},\bar{b})$	$\displaystyle=\sum_{\mathbf{G}}P_{\mathbf{k}}(\mathbf{r},\mathbf{G})C_{\mathbf{k}}(\mathbf{G},\bar{b})$		(20)

where the bar indicates that we are processing only a subset of the total number of bands. The phase matrix requires $\mathcal{O}[N]$ operations from the plane waves, regardless of the number of processors included. Projecting the wave functions is $\mathcal{O}[N^{2}]$ from the plane waves and bands, but the bands are distributed by processor. The summation over k-points is not counted in the estimation of computational cost because it decreases with volume and is usually 8 or 1. For a system with more than one site of interest, e.g., a disordered, liquid, or amorphous system, the number of sites, and therefore the number of local real-space grids, increases with volume as well. This means that the actual costs increase to $N^{2}$ and $N^{3}$.

To avoid $N^{3}$ scaling, we instead use a fast Fourier transform (FFT), followed by interpolation, and completed by applying the complex phase:

	$\displaystyle u_{\bar{b}\mathbf{k}}(\mathbf{x})$	$\displaystyle=\mathcal{FFT}\left[u_{\bar{b}\mathbf{k}}(\mathbf{G})\right]$
	$\displaystyle u_{\bar{b}\mathbf{k}}(\mathbf{r})$	$\displaystyle=\mathrm{Interp.}\left[u_{\bar{b}\mathbf{k}}(\mathbf{x})\right]$		(21)
	$\displaystyle\psi_{\bar{b}\mathbf{k}}(\mathbf{r})$	$\displaystyle=e^{i\mathbf{k}\cdot\mathbf{r}}u_{\bar{b}\mathbf{k}}(\mathbf{r}).$

The real-space grid $\mathbf{x}$ is defined as the Fourier transform dual of $\mathbf{G}$. We use 3rd-degree Lagrange polynomial interpolation to generate the wave functions on our desired points $\mathbf{r}$ from the FFT grid $\mathbf{x}$, aided by oversampling the FFT. The interpolants are cached to allow reuse within and between sites. The costs, including a factor of $N$ sites, are $\mathcal{O}[N^{2}\log N]$, $\mathcal{O}[N^{2}]$, and $\mathcal{O}[N^{2}]$, respectively. All three steps are independent over bands, k-points, and spins, providing good scaling with the number of processors.

To determine the break-even point between these two methods we must be more specific with the actual costs of each step. The $N^{3}$ term from method 1 is $N_{\mathbf{G}}N_{b}N_{\mathbf{r}}N_{i}$, where $i$ are the atomic sites. The $N^{2}\log N$ term from method 2 is $A_{F}N_{\mathbf{G}}N_{b}\log N_{\mathbf{G}}$, where $A_{F}$ is the FFT prefactor. Therefore, method 2 is preferable if $A_{F}\log N_{\mathbf{G}}<N_{\mathbf{r}}N_{i}$. Under the assumption that the logarithm of even a large number is about 10, method 2 is likely preferable, even for single-site calculations [60]. More sophisticated methods for Fourier interpolation onto irregular grids have been proposed in literature, e.g., [61] and references therein, but are not explored here.

The next step is augmenting the wave functions to recreate the all-electron character,

	$\displaystyle\psi^{\textit{ae}}_{\bar{b}\bar{\mathbf{k}}}(\mathbf{r})$	$\displaystyle=\psi^{\textit{ps}}_{\bar{b}\bar{\mathbf{k}}}(\mathbf{r})$		(22)
		$\displaystyle+\sum_{\nu,l,m}Y_{lm}(\hat{r})\left[\phi^{ae}_{\nu l}(\mathbf{r})-\phi^{ps}_{\nu l}(\mathbf{r})\right]\langle\phi^{ps}_{\nu l}Y_{lm}|\psi^{\textit{ps}}_{\bar{b}\bar{\mathbf{k}}}\rangle\;,$

where $\phi$ are the OPFs and $Y_{lm}$ are spherical harmonics. The local basis is substantially coarser than that used in the construction of the OPFs. Therefore we enforce unitarity by constructing the overlap matrix $A$,

$$A_{\nu\nu^{\prime};l}=\int_{0}^{r_{a}}drr^{2}\phi^{\textit{ps}}_{\nu l}(r)\phi^{\textit{ps}}_{\nu^{\prime}l}(r)\;,$$

(23)

where any deviation from the identity matrix is due to errors from using a coarser grid. The augmentation of Eq. 22 is modified,

	$\displaystyle\psi^{\textit{ae}}_{\bar{b}\bar{\mathbf{k}}}(\mathbf{r})=\psi^{\textit{ps}}_{\bar{b}\bar{\mathbf{k}}}(\mathbf{r})+\sum_{lm}$	$\displaystyle\sum_{\nu,\nu^{\prime}}Y_{lm}(\hat{r})\left[\phi^{ae}_{\nu l}(\mathbf{r})-\phi^{ps}_{\nu l}(\mathbf{r})\right]\times\quad$
	$\displaystyle\times$	$\displaystyle A^{-1}_{\nu\nu^{\prime};l}\langle\phi^{ps}_{\nu^{\prime}l}Y_{lm}|\psi^{\textit{ps}}_{\bar{b}\bar{\mathbf{k}}}\rangle,$		(24)

preserving unitarity.

Each process stores a copy of the OPFs and carries out the augmentation for its subset of bands and k-points. The scaling of this section goes as $\mathcal{O}[N^{2}]$ with a factors of $N_{i}$ sites and $N_{b}$ bands.

The RPA polarizability from Eq. 5 can be transformed from the two-time form to a convolution over energy, which can be carried out along the imaginary axis,

	$\displaystyle\chi_{0}$	$\displaystyle(\mathbf{r},\mathbf{r}^{\prime},\omega=0)=-i\sum_{\sigma}\int_{-\infty}^{\infty}\frac{dE}{2\pi}\,g_{\sigma}(\mathbf{r},\mathbf{r}^{\prime},E)g_{\sigma}(\mathbf{r}^{\prime},\mathbf{r},E)$
		$\displaystyle=\sum_{\sigma}\int_{-\infty}^{\infty}\frac{dt}{2\pi}\,g_{\sigma}(\mathbf{r},\mathbf{r}^{\prime},\mu+it)g_{\sigma}(\mathbf{r}^{\prime},\mathbf{r},\mu+it)$		(25)

where $\mu$ is chosen to be in the middle of the gap for insulators or at the Fermi level in metals to avoid poles in $g$. [With minimal approximation, a small energy $\Delta$ can be added in quadrature in metals to the difference between $\mu$ and the Kohn-Sham eigenvalue $\varepsilon_{b\mathbf{k}}$ according to $(\mu-\varepsilon_{b\mathbf{k}})\rightarrow\pm\sqrt{(\mu-\varepsilon_{b\mathbf{k}})^{2}+\Delta^{2}}$.] This same approach can be used for calculating the dynamic polarizability, $\omega\neq 0$. However, additional care is needed to avoid the poles in Green’s function $g(\omega+\mu+it)$ as $t\rightarrow 0$. The integral is replaced by a sum over an energy grid as outlined in appendix D.1.2.

In principle, partial Green’s functions could be constructed using the band and k-point distribution from the previous step. However, it is more efficient to redistribute the wave functions into blocks by $\mathbf{r}$ and site. The processors are split into groups such that each group works on its own site or set of sites. Within each group, the processors divide up the $\mathbf{r}$-points. This means that the wave functions are now distributed as $\psi_{b\mathbf{k}\sigma}(\bar{\mathbf{r}};\bar{\tau})$.

$$g^{(\tau)}_{\sigma}(\bar{\mathbf{r}},\mathbf{\bar{r}}^{\prime},\mu+it)=N_{k}^{-1}\sum_{\mathbf{k}}\sum_{b=1}^{N_{b}}\frac{\psi^{\dagger}_{b\mathbf{k}\sigma}(\bar{\mathbf{r}})\psi_{b\mathbf{k}\sigma}(\bar{\mathbf{r}}^{\prime})}{\mu+it-\varepsilon_{b\mathbf{k}\sigma}}$$

(26)

where $\mathbf{r}$ implicitly includes only the points for site $\tau$. If background communications are enabled, the majority of this data transfer takes place while the conversion process is on-going. Within a group of processors cooperatively working on a site $g^{(\tau)}$, the wave functions are shared. The scaling of this section goes as $\mathcal{O}[N^{2}]$ with a factors of $N_{i}$ sites and $N_{b}$ bands.

An important consideration in efficiently calculating the Green’s function is that it involves an outer-product of the wave functions. For each band and k-point, $N_{r}$ inputs are turned into $N_{r}^{2}$ outputs, with $2N_{r}^{2}$ operations. In a typical, small calculation the real-space grid has 1600 points per site, which means that at each frequency the Green’s function is just under 40 MB in size, making it too large to fit in the local cache of a typical modern processor. A naïve implementation would be limited by memory bandwidth. Instead, the wave functions are broken up, and the Green’s function is calculated by tiles in real space.

In the previous step, we calculated the irreducible polarizability $\chi^{0}$. The reducible polarizability is given by

$\displaystyle\chi=\left[1-\chi^{0}v\right]^{-1}\chi^{0},$

(27)

where $v$ is the Coulomb potential operator. We do this by projecting into a spherical basis

	$\displaystyle\chi^{0}(\mathbf{r},\mathbf{r}^{\prime})$	$\displaystyle=\sum_{lm}\sum_{l^{\prime}m^{\prime}}\chi^{0}_{lm;l^{\prime}m^{\prime}}(r,r^{\prime})Y_{lm}(\hat{r})Y^{*}_{l^{\prime}m^{\prime}}(\hat{r}^{\prime})$		(28)
	$\displaystyle v(\mathbf{r},\mathbf{r}^{\prime})$	$\displaystyle=\sum_{lm}\frac{4\pi}{2l+1}\frac{r_{<}^{l}}{r_{>}^{l+1}}Y_{lm}(\hat{r})Y^{*}_{lm}(\hat{r}^{\prime})$		(29)

where the Coulomb operator is diagonal in $l,m$. We can define

	$\displaystyle S_{lm;l^{\prime}m^{\prime}}(r,r^{\prime})$	$\displaystyle=\delta_{l,l^{\prime}}\delta_{m,m^{\prime}}\delta(r-r^{\prime})$		(30)
		$\displaystyle-\int\!dxx^{2}\,\chi^{0}_{lm;l^{\prime}m^{\prime}}(r,x)\frac{4\pi}{2l^{\prime}+1}\left[\frac{r_{<}^{l^{\prime}}}{r_{>}^{l^{\prime}+1}}\right]_{x,r^{\prime}}$

by taking advantage of the diagonal nature of $v$ in this basis. We therefore have

$\displaystyle\chi_{lm;l^{\prime}m^{\prime}}(r,r^{\prime})=S^{-1}_{lm;l^{\prime\prime}m^{\prime\prime}}(r,r^{\prime\prime})\chi^{0}_{l^{\prime\prime}m^{\prime\prime};l^{\prime}m^{\prime}}(r^{\prime\prime},r^{\prime}),$

(31)

where the matrices have dimension $N_{r}(2N_{l}+1)$.

In this basis the induced change in electron density from the short-range part of the core-hole potential $v^{(1)}$ is

	$\displaystyle\rho^{\textrm{ind}}(\mathbf{r})$	$\displaystyle=\sum_{lm}\rho_{lm}^{\textrm{ind}}(r)Y_{lm}(\hat{r})$		(32)
	$\displaystyle\rho_{lm}^{\textrm{ind}}(r)$	$\displaystyle=\sum_{l^{\prime}m^{\prime}}\int\!d^{3}r^{\prime}\chi_{lm;l^{\prime}m^{\prime}}(r,r^{\prime})v^{(1)}(r^{\prime})Y_{l^{\prime}m^{\prime}}(\hat{r}^{\prime})$
		$\displaystyle=\int dr^{\prime}r^{\prime 2}\chi_{lm;00}(r,r^{\prime})v^{(1)}(r^{\prime})$		(33)

The perturbing (core-hole) potential is taken to be spherical and therefore only the $l^{\prime}=0$ part of $\chi$ contributes. Giving a final, induced potential

	$\displaystyle v^{\textrm{ind}}(\mathbf{r})$	$\displaystyle=\sum_{lm}v_{lm}^{\textrm{ind}}(r)Y_{lm}(\hat{r})$
	$\displaystyle v_{lm}^{\textrm{ind}}(r)$	$\displaystyle=\int\!dr^{\prime}r^{\prime 2}\rho_{lm}^{\textrm{ind}}(r^{\prime})\frac{4\pi}{2l+1}\frac{r_{<}^{l}}{r_{>}^{l+1}}.$		(34)

By default, only the $l=0$ part of the response is calculated, according to

	$\displaystyle\bar{\chi}(r,r^{\prime})$	$\displaystyle=S^{-1}_{00;00}(r,r^{\prime\prime})\chi^{0}_{00;00}(r^{\prime\prime},r^{\prime})$
	$\displaystyle\bar{\rho}^{\textrm{ind}}(r)$	$\displaystyle=\int dr^{\prime}r^{\prime 2}\bar{\chi}(r,r^{\prime})v^{(1)}(r^{\prime})\,.$		(35)

The resulting induced potential is approximately the same as the $l=0$ component of the full induced density $\bar{\rho}\approx\rho_{00}$. For the core-hole potential, the strong localization means that $l=0$ component of the induced potential is dominant, and this approximation is reasonable. For valence calculations, including $l\leq 2$ was found to be sufficient. Because the dimension of $\chi$ is independent of system size, this section scales only with the number of sites $N_{i}$, linearly with system size or $\mathcal{O}[N]$.

For large cells, only a single k-point is required, and the electron orbitals can be calculated at the $\Gamma$-point. For systems with time-reversal symmetry the Bloch functions can be treated as real (instead of complex). This results a reduction of the required storage by half and substantial time savings in the DFT stage. A smaller reduction in runtime is also realized in the screening calculation as shown below in Table 1.

The calculation of the screening as outlined here is dominated by three steps: calculating the wave functions, projecting them onto the radial grid, and constructing the Green’s function and $\chi^{0}$. The first, calculating the electron orbitals using DFT, is carried out using the Quantum ESPRESSO code [22, 23, *espresso0]. We report the timing of the DFT step for completeness, however, we are focused on the two steps that are specific to the screening calculation.

To investigate the timing and scaling of screening calculations within ocean we use LiF (physical details and convergence parameters are given in appendix A). There are two classes of scaling that we are interested in. First there is system scaling, by which we mean the increase in run time with an increase in the system size. This highlights the inherent simulation size limits of our approach. We will also consider strong scaling, the change in execution time due to changing the number of processors. We have implemented two levels of parallelism for the screening calculations: internode mpi and shared memory openmp. The testbed for these calculations is a small cluster with 12 nodes. Each node has a dual-socket with 8 processors per CPU (16 per node) [62]. Each timing run was repeated 8 times, and the average value is reported (DFT calculations were run only once).

For the tests, we consider various supercells of LiF, from the unit cell $N_{F}=1$, to an $6\times 8\times 8$ supercell $N_{F}=384$, covering cell volumes from 110 a.u.³ to 42000 a.u.³ (6.2 nm³). For these runs, 32 bands per unit cell were included (12288 for $N_{F}=384$). For each supercell the screened Coulomb potential is calculated for all the fluorine sites. Each local real-space grid had 2624 points, and the Green’s functions were evaluated at 16 imaginary frequencies.

First we look at the scaling with system size or weak scaling. In Fig. 12 we show the run time of the projection ($\psi$) and construction of the Green’s function and polarizability ($g$ & $\chi$) steps as a function of super cell size. The straight lines on the log-log plot are $\alpha N^{2}$, where $\alpha$ is set by the timing of $N_{F}=125$. In this set of runs only $\Gamma$-point sampling of the Brillouin zone was used. The $\approx\mathcal{O}[N^{2}]$ growth in calculating the screening is evident by the linear plots, though overhead or inefficiencies dominate the timing of the smallest run. Additionally, for large systems the ‘$g$ & $\chi$’ diverges slightly from the expected $N^{2}$ behavior which may be indicative of poor cache reuse or other bottlenecks for large system sizes. We also include the timing of the DFT portion which follows $\alpha N^{3}$. 128 processors across 8 nodes were used for each run, using 128 mpi processes.

The timing information for the range of systems from $N_{F}=1$ to $N_{F}=384$ is shown in Table 1. To give a full picture of the scaling, three different settings for the k-point sampling were used: finite sampling on a $2^{3}$ k-point mesh, single k-point sampling, and $\Gamma$-point sampling. Above $N_{F}=27$, a single k-point is sufficient to sample the Brillouin zone and gives the same results as the $2^{3}$ sampling (but more quickly). The purpose of timing single k-point runs in addition to the $\Gamma$-point runs is to distinguish the changes due to the reduction in k-point sampling from 8 to 1 from the changes in moving between complex and real Bloch functions. The single k-point is taken at $(\nicefrac{{1}}{{8}},\nicefrac{{2}}{{8}},\nicefrac{{3}}{{8}})$.

Next, we present the change in run time with changing processor number or strong scaling, Fig. 13. Ideally, doubling the number of processors used in a calculation will halve the runtime. Longer than expected runtimes may result from serial sections of the code or communication overhead. Shorter than expected times may result from better data caching due to each processor working on a smaller data set. Here we plot the data as the efficiency $E$ as a function of the number of processors $N$,

$$E(N)=\frac{N_{0}}{N}\,\frac{t(N_{0})}{t(N)}\times 100\%$$

(36)

where the efficiency is normalized to the run time with $N_{0}$ processors. Ideal scaling is given by an efficiency of 100 %. The efficiency is the measure of merit for planning high-throughput calculations. In high-throughput calculations the available hardware resources can be divided between many different calculations, and the runtime of any single calculation should be balanced against the runtime of the complete dataset.

In Fig. 13(a) we show that for a moderately sized system $N_{F}=64$, there is a drop-off in efficiency above 16 processors. In part this is a reflection of the structure of our computer cluster, where each addition of 16 processors increases the number of nodes in the calculation by 1. While the efficiency is quite poor running on 160 processors, the runtime is also very brief. The average time is 41.1 sec. compared to an idea time of 27.3 sec. The larger systems show better scaling.

Lastly, we examine the omp parallelism by repeating the $\Gamma$-centered calculations this time including thread-level parallelism. The total number of processors is held fixed at 128, but divided between mpi and omp with 1, 2, 4, and 8 omp threads per mpi process. The results are shown in Table 2 and for three of system sizes in Fig. 13(b). We find relatively uniform performance across the first three processor arrangements, but a drop-off in performance using 8 omp threads. This drop-off indicates an opportunity for further code refinement to better support higher-levels of openmp parallelism. In Table 2 it can be seen that this inefficiency is primarily the result of poor scaling of the ‘$\psi$’ section.

In the current implementation there is no re-use of the Green’s functions between sites. The use of site-centered radial grids makes it unlikely that a given $(\mathbf{r,r}^{\prime})$ pair of one site will exist in the grid of another. However, it reasonable to expect that many point pairs will be nearly shared, e.g., for sites $\alpha$ and $\beta$ that $\mathbf{r}_{\alpha}\approx\mathbf{r}_{\beta}$ and $\mathbf{r}^{\prime}_{\alpha}\approx\mathbf{r}^{\prime}_{\beta}$. In the future, the construction of the grids can be relaxed to maximize the overlaps, decreasing the computational cost of generating the Green’s functions. This would be especially helpful in the case of valence calculations where the site density is high.

Future improvements to the scalability with system size must focus on generating the electron wave functions. For medium to large system sizes, most of the time is in calculating the screening is spent in the DFT (see Sec. VI.7). This is exacerbated by the need for unoccupied states in the calculating the Green’s function. Several methods have been proposed to reduce the number of unoccupied states.

One option is to directly replace part of the sum over unoccupied states. The effective-energy techniques replace the energy denominator in the sum over unoccupied states [63, 64, 65]. The completeness of the Bloch functions then allows the sum over unoccupied bands to be replaced with the identity minus a sum over occupied bands. However, these approaches differ in two main ways from ours, not including our local approximation. First, the energy convolution is carried out analytically, and the RPA polarizability is constructed via sums over states. In our approach the Green’s functions are built via a sum over states and the convolution is carried out numerically. Second, the effective-energy approaches are formulated in reciprocal space, which has the advantage of a straightforward approximation for the effective energy. An easier approach might be to approximate the neglected high-energy bands as plane waves [66, 67, 68, 69].

Alternatively, the induced density response and therefore the screened potential can be more directly calculated using the linear-response Sternheimer equation approach [6, 70] or eigenvalue decomposition of the polarizability matrix [71, 7]. While these approaches only require the occupied orbitals, they maintain an unfavorable $N^{4}$ scaling with system size. Better scaling might be achievable by adapting these approaches to determine only the local response.

Lastly, the model dielectric function used to screen $v_{2}$ (Eq. II.2) was designed for bulk systems. In the case of highly anisotropic systems, like a surface or interface, this may result in slower convergence with respect to the shell radius $R_{S}$, requiring a larger real-space RPA calculation. This could be addressed by modifying the screening model or alleviated through a judicious choice of model parameters, namely choosing the average density $\rho_{0}$ to more accurately reflect the bulk material (see model details in Appendix B).

The screened Coulomb interaction $W$ has many uses in calculations of condensed matter systems other than the use presented here of the direct interaction in the the BSE. One such application is in self-energy calculations using the GW method which requires evaluating the self-energy operator [72],

$$\Sigma(\mathbf{r},\mathbf{r}^{\prime},E)=i\int\!\frac{d\omega}{2\pi}\;G(\mathbf{r},\mathbf{r}^{\prime},E-\omega)W(\mathbf{r},\mathbf{r}^{\prime},\omega)\quad.$$

(37)

The local screening approach outlined here can be used to efficiently generate $W$ with RPA quality for small distances $|\mathbf{r}-\mathbf{r}^{\prime}|$, just as was shown for valence BSE calculations. The terms in the GW calculations are significantly less localized than core-level excitations, making discrepancies in the short-range part of $W$ less important, and augmentation of the orbitals may not be necessary. However, in transition metals the d orbitals often drive important characteristics of the electronic behavior, forming the top of the valence bands, the bottom of the conduction bands, or both. From atomic calculations, it can be seen that the d orbitals overlap significantly with the semi-core orbitals of the same principle quantum number, e.g., the 3d with the 3s and 3p. High-accuracy calculations involving localized d orbitals may require accurately correcting the nodal structure of the s and p wave functions in much the same manner as the we have shown for core-level spectroscopy. Several GW studies have pointed out discrepancies from using pseudopotentials [39] or excluding semi-core states [73, 74]. In the present approach only static screening was implemented. However, the contour integral in Eq. 25 can be modified to calculate $\omega\neq 0$, and the fundamental scaling of frequency-dependent screening remains the same as that of the static case. In addition to standard valence- and conduction-band self-energy calculations, an adaptation of this method could be applicable for determining accurate core-level binding energies [75].

The local, real-space screening might also be useful for phonon calculations. Within the harmonic approximation, the phonons of a system can be fully described by the dynamical matrix. The elements of the dynamical matrix are proportional to the derivative of the force on atom $a$ with respect to changes in the position of $a^{\prime}$. This is equivalent to the second derivative of the total energy with respect to the displacement of both. The elements of the dynamical matrix can be calculated using density-functional perturbation theory [76, 77, 78]. Alternatively, the dielectric response can be used since the polarizability describes how the electron density will change in response to a change in the potential, in this case the motion of the atomic nuclei. Care would be required as the polarizability gives the density change in response to a local perturbing potential, but standard pseudopotentials include non-local terms [79].

The screening calculations in this paper have been carried out only within the RPA approximation. From a many-body perturbation theory perspective, the RPA is the lowest-order diagram for the polarizability. Higher-order approaches treat interactions between the electron and hole lines in the RPA, or, equivalently, add a vertex correction. Unfortunately, additional interaction or vertex terms increase the computational cost and scale worse with system size. Our local, real-space approach is an ideal testbed for investigating higher-order approaches because the increase in scaling complexity applies to the dimension of our local dielectric response function which is small and independent of the system size.

As an example, we have implemented the vertex correction given by the adiabatic local-density approximation (ALDA). As shown by Del Sole, Reining, and Godby [80], if the first-order approximation to the one-electron self-energy is taken to be the local exchange-correlation potential

$\displaystyle\Sigma(1,2)=\delta(1,2)v_{\textit{xc}}(1)$

(38)

then the reducible polarizability (and screened interaction $W$) undergo a relatively simple transformation. Repeating Eq. II.1,

	$\displaystyle\chi$	$\displaystyle=\left(1-\chi_{0}v\right)^{-1}\chi_{0}$
	$\displaystyle\tilde{\chi}$	$\displaystyle=\left(1-\chi_{0}(v+f_{\textit{xc}})\right)^{-1}\chi_{0}$		(39)

where $f_{\textit{xc}}$ is the derivative of the exchange-correlation potential with respect to the density, $f_{\textit{xc}}(1,2)=\delta v_{\textit{xc}}(1)/\delta n(2)$, and $\tilde{\chi}$ is the ALDA polarizability. The use of an ALDA kernel has been investigated within the GW approximation [80, 81, 82] and for valence BSE calculations of small molecules [83].

Within the ALDA, $f_{\textit{xc}}$ is a contact interaction, and the expression in Eq. 39 is easily evaluated using the ocean code as outlined in section VI.5. In the real-space basis $f_{\textit{xc}}$ is diagonal and can be written

$$f_{\textit{xc}}(\mathbf{r},\mathbf{r}^{\prime})=\delta(\mathbf{r-r}^{\prime})\left.\frac{dv_{\textit{xc}}(n)}{dn}\right|_{n=n(\mathbf{r})}$$

(40)

where $v_{\textit{xc}}$ is the LDA exchange-correlation potential and is evaluated at the density $n$ at position $\mathbf{r}$. The electron density $n(\mathbf{r})$ is taken from the initial DFT calculation used to generate the electron orbitals for the screening. We start with the Perdew-Zunger parameterization for the exchange [84] and Vosko, Wilk and Nusair parameterization of the correlation energy [85] within the local-density approximation fit to the data of Ceperley and Alder [86]. We calculate $f_{xc}$ directly as the second derivative of the exchange-correlation energy with respect to the density using a 5-point finite difference using density differences of 0.01 e^- per a.u.³. Spin-polarized calculations beyond the RPA are not yet supported, but can be included using this same scheme.

Once again looking at h-BN, we can examine how calculating the polarizability with the ALDA instead of RPA changes the XAS. In general, the ALDA results in a stronger induced potential [shown for LiF in Fig. 16(a)]. This in turn leads to a weaker core-hole potential and correspondingly weaker excitonic effects. In Fig. 14 we show the nitrogen K-edge XAS using both the RPA and ALDA for the screening. The BSE spectrum calculated using the ALDA is substantially different from the RPA result, but only in the near-edge region, within about 15 eV of the onset. The small differences at higher energies would be hidden by broadening if the calculation included the effects of the electron self-energy and vibrational disorder. While electron energy loss spectroscopy (EELS) taken within the dipole limit should probe the same excitations as XAS, it is seen in h-BN that there are large discrepancies between the two [43]. In part, this is due to the different systematic errors such as surface sensitivity and self-absorption effects affecting XAS versus EELS. The ALDA screening appears superior to the RPA when comparing to the EELS data, but the RPA appears superior when comparing to the XAS measurement. A broad survey of materials and careful quantification of experimental uncertainties is necessary to establish the general applicability of one approximation or the other and should be the subject of future work.

While the only vertex correction that has been implemented in ocean is the ALDA, an extension to semi-local exchange-correlation kernels is straightforward, requiring only the additional knowledge of the density gradients. Because $f_{\textit{xc}}$ is formed explicitly in real-space, the formulation of Eq. 39 is also compatible with non-local exchange-correlation potentials. This would require construction of $f_{\textit{xc}}$ as a real-space matrix instead of the diagonal form (Eq. 40). However, the response is still localized and in response to a local perturbing potential with the long-range response handled by a model. Therefore, any non-local $f_{\textit{xc}}$ must also be of limited range.