Stochastic density functional theory (sDFT) is becoming a valuable tool for studying groundstate properties of extended materials. The computational complexity of describing the Kohn–Sham orbitals is replaced by introducing a set of random (stochastic) orbitals leading to linear and often sublinear scaling of certain groundstate observables at the account of introducing a statistical error. Schemes to reduce the noise are essential, for example, for determining the structure using the forces obtained from sDFT. Recently, we have introduced two embedding schemes to mitigate the statistical fluctuations in the electron density and resultant forces on the nuclei. Both techniques were based on fragmenting the system either in real space or slicing the occupied space into energy windows, allowing for a significant reduction in the statistical fluctuations. For chemical accuracy, further reduction of the noise is required, which could be achieved by increasing the number of stochastic orbitals. However, the convergence is relatively slow as the statistical error scales as $1/ N \chi $ according to the central limit theorem, where N_{χ} is the number of random orbitals. In this paper, we combined the embedding schemes mentioned above and introduced a new approach that builds on overlapped fragments and energy windows. The new approach significantly lowers the noise for groundstate properties, such as the electron density, total energy, and forces on the nuclei, as demonstrated for a Gcenter in bulk silicon.
I. INTRODUCTION
Kohn–Sham (KS) density functional theory^{1,2} (DFT) is widely used to study a wide range of systems due to its capability of quantitatively predicting groundstate properties at a moderate computational cost of $O ( N e 3 ) $, where N_{e} is the number of electrons. While this moderate scaling allows for an efficient description of the ground state of molecules and bulk structures with periodic boundary conditions, the application to systems containing 10^{4} electrons or more, such as nanostructures,^{3} complex materials,^{4} and large biomolecules,^{5} is still a severe challenge for today’s DFT implementations. Linearscaling methods for DFT based on dividing the entire system into subsystems^{6–8} require proper treatment of the boundaries between the fragments.^{8–13} Another approach to realize linear scaling relies on the “locality” of the density matrix (DM).^{14–20} These approaches have received growing interest in recent years.^{21–25}
We have recently introduced an alternative linearscaling approach to DFT, which does not rely on the partitioning of the system into subsystems nor depend on the sparsity of the density matrix.^{26} Instead, it utilizes stochastic orbitals, which are random linear combinations of deterministic KS orbitals, to calculate the electron density and other groundstate properties. In practice, the required number of stochastic orbitals does not increase with the system size for evaluating many groundstate properties,^{26,27} leading to linear scaling or even sublinear scaling DFT.^{26,27} In stochastic DFT, linear scaling is achieved at the cost of introducing a statistical error in the density and related observables, which, according to the central limit theorem, decrease rather slowly with the number of stochastic orbitals, N_{χ}, limiting the efficiency and accuracy of the method. Therefore, developing noise reduction schemes for sDFT is essential for achieving chemical accuracy without the need to dramatically increase N_{χ}.
One approach for reducing the noise in sDFT is based on dividing the entire system into fragments. The entire system’s density is then given as a sum of the fragment densities and a correction term sampled using stochastic orbitals. When the sum of the fragment densities provides a good approximation of the total system’s density, the correction term is small, leading to significant reductions in the noise for the electron density, energy, and forces on the nuclei. This approach has been illustrated for systems with open boundary conditions^{28–30} as well as for periodic boundary conditions.^{27,31} For the latter case, we used overlapped fragments to ensure a reasonable estimate of both the density and the density matrix [this approach was referred to as “overlapped embeddedfragmented stochastic density functional theory” (oefsDFT).^{27}] Recently, we introduced an alternative technique to mitigate the statistical noise, referred to as “energywindow sDFT” (ewsDFT),^{32} where the occupied space was divided into energyresolved subspaces (“energy windows”) and the contribution to the density for each window can be calculated simultaneously. This method reduces the statistical noise in the density and the nuclei forces, but not in the total electronic energy.^{32}
In this paper, we combine the overlapped embeddedfragmented scheme with the energy window scheme. Noise reduction is obtained by projecting both the system density matrix and the fragment density matrix onto fixed energy windows. The proposed energy window embeddedfragmented stochastic DFT (ewefsDFT) approach reduces the noise in the electron density, total energy, and forces on the nuclei and the total computational time by more than an order of magnitude compared to ewsDFT or oefsDFT, as illustrated for a Gcenter embedded in bulk silicon. The noise reduction is crucial for obtaining structural information with chemical accuracy using only several tens of stochastic orbitals, as will be shown in a proceeding publication.^{33}
This manuscript is organized as follows: In Sec. II, we briefly review the sDFT. In Sec. III, we present the oefsDFT and ewsDFT methods, both central to the development of the current noise reduction scheme. In Sec. IV, we provide the details of the proposed ewefsDFT and a summary of the algorithm. Assessment of the new approach for a challenging Gcenter embedded in bulk silicon is presented in Sec. V alongside a discussion of the computational complexity and cost of the ewefsDFT. Finally, in Sec. VI, we summarize the main developments.
II. STOCHASTIC DENSITY FUNCTIONAL THEORY
Consider an extended system described by KSDFT, with a KS Hamiltonian ($ h \u0302 K S $) given by
where $ t \u0302 $, $ v \u0302 nl $, $ v \u0302 loc $, $ v \u0302 H [ \rho ] $, and $ v \u0302 xc [ \rho ] $ are the operators of the kinetic energy, the nonlocal pseudopotential energy, the local pseudopotential energy, the Hartree energy, and the exchange–correlation energy, respectively. The Hartree and exchange–correlation terms depend on the electron density, $\rho r $, which is formally given by (we assume closedshell and ignore spin–orbit couplings for simplicity)
where $ \rho \u0302 = \theta \beta ( h \u0302 KS , \mu ) $ is the onebody density matrix and θ_{β}(x, μ) = 1/(1 + e^{β(x−μ)}) is the Fermi–Dirac distribution function parameterized by the inverse temperature (β) and the chemical potential (μ) tuned to give the number of electrons, N_{e} = ∫drρ(r). Other smooth functions to approximate a step function can also be used instead of θ_{β}(x, μ). In KSDFT, the electron density can also be written in terms of the KS orbitals (eigenstates of the KS Hamiltonian), ϕ_{i}(r),
where N_{occ} is the number of occupied orbitals.
In sDFT, the trace in Eq. (2) is replaced by averaging the expectation value of $ \theta \beta ( h \u0302 KS , \mu ) \delta ( r \u2212 r \u0302 ) $,
where χ⟩ is a stochastic orbital and ⟨⋯⟩_{χ} implies averaging over an ensemble of stochastic orbitals. The stochastic orbitals are represented on a realspace grid with N_{g} grid points; each grid point is assigned a random value $\xb11/ \Delta V $, where ΔV = V/N_{g} is the volume element and V is the volume of the supercell. Equation (4) can be rewritten in a compact form as
where ξ⟩ is a projected stochastic orbital, which is given as
The projection of the stochastic orbitals onto the occupied space is obtained by expanding $ \theta \beta ( h \u0302 KS , \mu ) $ in Chebyshev polynomials,^{34,35}
where N_{c} is the length of the Chebyshev polynomial expansion, a_{n}(μ, β) are the expansion coefficients, and T_{n} are the Chebyshev polynomials of order n.
Groundstate observables corresponding to any onebody operator, $ O \u0302 $, can be evaluated using a similar stochastic trace formula,
Since the exact electron density can only be recovered by averaging infinitely many stochastic orbitals, estimates for O result in statistical errors that decrease as $ N \chi \u2212 1 / 2 $ according to the central limit theorem, where, as before, N_{χ} is the number of stochastic orbitals. To achieve chemical accuracy for the electron density and the forces on the nuclei, N_{χ} may need to exceed 10^{3} orbitals, limiting the efficiency of sDFT. The need to develop noise reduction schemes is clear and would extend the range of system sizes that can be studied routinely using sDFT.
III. NOISE REDUCTION SCHEMES IN STOCHASTIC DFT
A. Overlapped embeddedfragmented stochastic DFT
Significant reduction in the statistical error can be achieved by introducing a reference system that provides a reasonable approximation to the electron density and can be calculated within KSDFT. The total electron density is then given as a sum of the reference system electron density and a small correction term obtained stochastically.^{27–30} In this section, we will briefly review the most recent developments based on an overlapped embeddedfragmented stochastic DFT (oefsDFT), which is central to the proposed ewefsDFT. Full details of the approach can be found elsewhere.^{27}
In oefsDFT, the supercell is divided into fragments referred to as “core regions” (see Fig. 1 for an illustration) C_{f} wrapped by “buffer regions” (B_{f}) to form dressed fragments (D_{f} = C_{f} ∪ B_{f}), where f is the fragment index. The fragment density matrix, $ \rho \u0302 f $, is given by
for r ∈ C_{f}. In the above equation, φ^{f}(r) are the KS orbitals for fragment f obtained from a deterministic KSDFT approach and $ N occ f $ is the total number of occupied orbitals for the fth dressed fragment. Using the above equation, the total electron density can be evaluated as follows:
where the fragment electron density is $ \rho f ( r ) = \u2211 i = 1 N occ f  \phi i f ( r )  2 $, $ \xi f ( r ) = \u2211 i = 1 N occ f \phi i f ( r ) \u27e8 \phi i f  \chi \u27e9 D f $, and $ \u27e8 \phi i f  \chi \u27e9 D f = \u222b D f dr \phi i f ( r ) * \chi ( r ) $. We use the relationship $ \rho f ( r ) = \u27e8 r  \rho \u0302 f \rho \u0302 f \u22a4  r \u27e9 $ in Eq. (10) since KSDFT methods in the 0 K limit are adopted to calculate fragment KS orbitals.
In the limit N_{χ} → ∞, the first and last terms on the right hand side of Eq. (10) cancel, while the remaining term converges to the deterministic electron density. For a finite set of stochastic orbitals, the noise in the second term on the right hand side of Eq. (10) roughly cancels that in the last term, as long as the reference system density matrix provides a reasonable approximation to that of the full system, thereby leading to a significant reduction in the statistical error.^{27,28} Finally, we would like to note that the fragmentation scheme described above is similar to another fragmentationbased DFT approach.^{8,36} However, unlike divideandconquer based methods, the stochastic approach does not introduce systematic errors resulting from the fragmentation in the limit N_{χ} → ∞.
B. Energy window stochastic DFT
A reduction in the statistical fluctuations can also be achieved using another scheme, based on partitioning the occupied space into “energy windows.”^{32} In this approach, rather than projecting χ⟩ onto the occupied space using Eq. (6), we divide the occupied space into energy windows, and χ⟩ is projected onto each window using a set of projectors, $ P \u0302 1 ,\u2026, P \u0302 N w $. Here, the projector $ P \u0302 w $ is defined as
for 1 ≤ w ≤ N_{w}, where $\u2212\u221e= \epsilon 0 < \epsilon 1 <\cdots < \epsilon N w =\mu $. In ewsDFT, the electron density is given by the sum of all projected densities,
where $ \xi ( w ) \u3009= P \u0302 w \chi \u3009$ is a projected stochastic orbital for window w. ξ^{(w)}⟩ are calculated simultaneously with a single Chebyshev expansion, i.e.,
where N_{c} is the length of the Chebyshev polynomial (chosen sufficiently large to converge the expansion), $ b n ( w ) ( \epsilon w , \epsilon w \u2212 1 ) $ is the expansion coefficient (depends on the window), and T_{n}(x) is the Chebyshev polynomial of order n. The variance of the electron density in this scheme is $8 \u2211 w = 1 N w \rho ( w ) ( r ) 2 $, which is smaller than the variance in sDFT given by $8 \u2211 w = 1 N w \rho ( w ) ( r ) 2 $.^{32}
IV. ENERGY WINDOW EMBEDDEDFRAGMENTED STOCHASTIC DFT
Combining the energy window approach with the fragmentation approach results in the following expression for the electron density at a grid point r:
In the above equation, the projection operators on the energy windows are defined as
where $ { \epsilon } \u2261 \epsilon 0 \u2026 \epsilon N w \u2212 1 $ (ɛ_{0} = −∞) define the boundaries of the energy windows. We would like to signify the differences between the projection operators used for ewsDFT [cf. Eq. (11)] and those used in the current ewefsDFT approach [cf. Eq. (15)]. In ewsDFT, $ \u2211 w = 1 N w P \u0302 w $ equals the density matrix $ \rho \u0302 $, while in ewefsDFT, it equals the unit operator, $ I \u0302 $. Furthermore, in ewsDFT, the highest energy window is set to $ \epsilon N w =\mu $, while in ewefsDFT, the energy windows are held fixed (but not necessarily equally spaced) for the entire selfconsistent procedure and are chosen to be independent of the chemical potential, μ. The use of a fixed set of energy windows simplifies the onthefly calculations of chemical potential [see Eq. (18)].
In Eq. (14), the action of $ \rho \u0302 P \u0302 w $ and $ P \u0302 w $ on χ⟩ is obtained using a proper Chebyshev series,
Finally, as before, the density of each fragment is given by $ \rho f ( r ) = \u2211 i = 1 N occ f  \phi i f ( r )  2 $ and the stochastic projected orbitals for each fragment are given by
For a fixed set of energy windows, ɛ_{w}, the chemical potential (μ) can be obtained by solving
where $ \u222b C f dr$ imply that the integrals are performed in real space in region r ∈ C_{f}. In the above equation, $ \u27e8 \chi  \rho \u0302 ( \mu )  \chi \u27e9 $ is evaluated by expanding $ \rho \u0302 $ in a Chebyshev series $ \u27e8 \chi  \rho \u0302 ( \mu )  \chi \u27e9 = \u2211 n = 0 N c c n ( \mu ) \u27e8 \chi  T n ( h \u0302 KS )  \chi \u27e9 $ and the chemical potential is determined by solving for N(μ^{*}) = N_{e}, where N_{e} is the total number of electrons in the system.
Similar to the electron density, other groundstate observables such as the kinetic energy,
or the nonlocal pseudopotential energy,
or the nonlocal pseudopotential contribution to the forces on the nuclei,
is expressed in ewofsDFT as a sum of three terms, where the first term and the last term cancel each other in the limit N_{χ} → ∞ (see the supplementary material).
The proposed ewefsDFT method to reduce the noise in the density, energy, and forces on the nuclei can be summarized as follows:

Generate the KS orbitals $ { \phi i f ( r ) } $ for each dressed fragment and $ \rho f ( r ) = \u2211 i = 1 N o c c f  \phi i f ( r )  2 $ using a deterministic DFT. ρ(r) = 2∑_{f}ρ_{f}(r) is used as the initial electron density guess.

For each stochastic orbital χ(r) (defined above), calculate and store on the grid the projected stochastic orbital $ \zeta ( w ) ( r ) = \u27e8 r  P \u0302 w  \chi \u27e9 $ and also store the Chebyshev moments, $ \u27e8 \chi  T n ( h \u0302 KS )  \chi \u27e9 $.

For each window and for each stochastic orbital, generate and store on the grid $ \xi f ( w ) ( r ) = \u2211 i = 1 N occ f \phi i f ( r ) \u27e8 \phi i f  \xi ( w ) \u27e9 D f $.

Solve for μ^{*} (N(μ^{*}) = N_{e}) with the regula falsi method using the Chebyshev moments ($ \u27e8 \chi  T n ( h \u0302 KS )  \chi \u27e9 $), $ \xi f ( w ) ( r ) $, and ρ_{f}(r).^{30}

For each window and for each stochastic orbital, generate and store on the grid the stochastic projected orbitals $ \zeta ( w ) ( r ) = \u27e8 r  \rho \u0302 ( \mu * ) P \u0302 w  \chi \u27e9 $ using the chemical potential determined in the previous step.

Generate and store the electron density ρ(r) using Eq. (14) with all stochastic orbitals.

Update the density and the KS Hamiltonian using the iterative subspace (DIIS) method^{37} and repeat the above steps until selfconsistency is achieved, using the same random number seed.
V. APPLICATION TO GCENTER IN BULK SILICON
We demonstrate the ewefsDFT for a low Gcenter defect concentration in bulk silicon.^{38,39} We focus on the Atype^{38} Gcenter impurity embedded within a Si_{512} supercell (see the structure in Fig. 2). We performed Γ point DFT calculations using the Perdew–Burke–Ernzerhof (PBE)^{40} functional with Troullier–Martins normconserving pseudopotentials^{41} in the Kleinman–Bylander form.^{42} As a result of localized orbitals around the carbon atoms, a 40 Ry wave function cutoff (80 Ry for the density cutoff) was used, corresponding to the realspace grid spacing of 0.18 Å. Ingap states require a large β ≈ 900 in inverse Hartree to sufficiently converge the groundstate properties with respect to the electron temperature. 80 stochastic orbitals were used in both ewefsDFT and oefsDFT. The dressed fragments were Si_{64} with periodic boundary conditions, while the size of each core region was Si_{8}. We used 41 energy windows in ewefsDFT with fixed window boundaries chosen such that each window has roughly the same number of KS orbitals, thereby lowering the statistical noise. This requires an estimate of the density of states, which can be obtained directly from the projection of the stochastic orbitals^{26} or approximated from the density of states of the fragments or the pristine structure. For each single point calculation, SCF convergence was achieved when the difference of energy per electron between two consecutive iterations is below 10^{−8} Hartree (see Fig. 3). The fragment density was used as an initial guess for the SCF calculations.
A. Results
In Fig. 4, we assess the accuracy of ewefsDFT for the electron density (upper panel) and the standard deviation in the electron density (lower panel) for selected positions in the vicinity of the Gcenter. The density and standard deviation were calculated from 5 independent ewefsDFT or oefsDFT runs, with 80 stochastic orbitals for each run. For clarity, we plot the absolute value of the electron density difference between the stochastic and the deterministic calculations. Both the deviations from the deterministic electron density and the standard deviation obtained by the ewefsDFT (shown in red) are significantly smaller than the corresponding oefsDFT results (shown in blue). The noise of the electron density is significantly smaller, by approximately, a factor of 5 when compared to oefsDFT.
In Table I, we list the kinetic energy (E_{k}), the nonlocal pseudopotential energy (E_{nl}), the Hartree energy (E_{H}), the local pseudopotential energy (E_{loc}), the exchange–correlation energy (E_{xc}), and the total energy (E_{tot}), all per electron. The reference deterministic calculation converges (on the grid) to all significant digits shown. In parentheses, we provide the standard error, which is significantly smaller in ewefsDFT compared to oefsDFT for all quantities. We find that the standard error in the total energy per electron decreased by a factor of ≈3 when 41 windows were used. The total energy per electron in both ewefsDFT and oefsDFT are slightly outside one standard deviation from the deterministic DFT result.^{30,31} We note in passing that the ewsDFT approach (without fragmentation) does not reduce the noise in the total energy per electron, as discussed previously.^{32}
Method .  E_{k}/N_{e} .  E_{nl}/N_{e} .  E_{H}/N_{e} .  E_{loc} .  E_{xc} .  E_{tot}/N_{e} . 

dDFT  10.1767  5.7871  2.0194  −8.7085  −8.1027  −26.8197 
oefsDFT  10.1735(35)  5.7882(36)  2.0222(19)  −8.7074(40)  −8.1052(8)  −26.8205(11) 
ewefsDFT  10.1778(11)  5.7874(4)  2.0193(6)  −8.7086(11)  −8.1028(3)  −26.8186(4) 
Method .  E_{k}/N_{e} .  E_{nl}/N_{e} .  E_{H}/N_{e} .  E_{loc} .  E_{xc} .  E_{tot}/N_{e} . 

dDFT  10.1767  5.7871  2.0194  −8.7085  −8.1027  −26.8197 
oefsDFT  10.1735(35)  5.7882(36)  2.0222(19)  −8.7074(40)  −8.1052(8)  −26.8205(11) 
ewefsDFT  10.1778(11)  5.7874(4)  2.0193(6)  −8.7086(11)  −8.1028(3)  −26.8186(4) 
In Fig. 5, we plot the force on the nuclei along the xdirection for selected atoms obtained by the ewefsDFT (upper panel) and the oefsDFT (lower panel). Error bars indicate the standard deviation for each force. Clearly, the statistical fluctuations are significantly smaller for ewefsDFT compared to oefsDFT. In order to estimate the overall noise reduction efficiency, we averaged the standard deviations of F_{x} (nuclei force along the x axis) over all atoms. The averaged standard deviation of F_{x} is ≈0.53 eV/Å for oefsDFT and ≈0.09 eV/Å for ewefsDFT. Similar results were also obtained for other force components F_{y} and F_{z}, which implies that ewefsDFT standard deviations of nuclei forces are about a factor of 6 smaller than those in oefsDFT. In other words, to achieve a similar noise level in oefsDFT would require ≈30 times more stochastic orbitals. No bias was observed for forces on the nuclei in ewefsDFT.
B. Computational cost
In Table II, we summarize the total wall time and the contribution from the Chebyshev moments, $ \u27e8 \chi  T n ( h \u0302 KS )  \chi \u27e9 $, and from the projections of χ⟩ onto the occupied space and the energy windows. The number of stochastic orbitals, the size of the core and dressed fragments, the grid size, and all other parameters were the same for both approaches. The overall wall time seems to be very similar, comparing ewefsDFT and oefsDFT methods. Each SCF iteration is ≈50% longer in ewefsDFT, but the number of SCF iterations required to achieve a similar convergence is smaller in ewefsDFT, resulting in similar wall times. Note that the statistical error in ewefsDFT is much smaller than oefsDFT. To achieve similar statistical errors in oefsDFT, it would result in wall times that are roughly 30 times longer than ewefsDFT.
.  oefsDFT .  ewefsDFT . 

Wall time (h)  25.11  25.22 
Number of SCF iterations  25.6  16.4 
Time for one SCF iteration (h)  0.98  1.54 
Time for calculating $ \u27e8 \chi  T n ( h \u0302 KS )  \chi \u27e9 $ (h)  0.32  0.73 
Time for generating ξ⟩ or ξ^{(w)}⟩ (h)  0.64  0.72 
Time for projecting ζ^{(w)}⟩ with $ \rho \u0302 frag $ (h)  N/A  0.04 
.  oefsDFT .  ewefsDFT . 

Wall time (h)  25.11  25.22 
Number of SCF iterations  25.6  16.4 
Time for one SCF iteration (h)  0.98  1.54 
Time for calculating $ \u27e8 \chi  T n ( h \u0302 KS )  \chi \u27e9 $ (h)  0.32  0.73 
Time for generating ξ⟩ or ξ^{(w)}⟩ (h)  0.64  0.72 
Time for projecting ζ^{(w)}⟩ with $ \rho \u0302 frag $ (h)  N/A  0.04 
The main difference between the two methods is the computational time for generating the Chebyshev moments, $ \u27e8 \chi  T n ( h \u0302 KS )  \chi \u27e9 $. In oefsDFT, we used the relation $ T 2 n =2 T n 2 \u22121$ to evaluate $ T n ( h \u0302 KS ) \chi \u3009$ for n > N_{c}/2, thereby computing only 1/2 the number of moments. This relation cannot be used in ewefsDFT since each stochastic orbital is projected onto all energy windows to generate ξ^{(w)}⟩ and ζ^{(w)}⟩. Other smaller contributions to the computational time difference between the two methods can be traced to the need to generate N_{w} projected stochastic orbitals in ewefsDFT compared to only one projected orbital in oefsDFT, resulting in ≈20% increase in generating ξ⟩ vs ξ^{(w)}⟩. In addition, in ewefsDFT, one has to compute $ \xi f ( w ) \u3009$ with the fragment density matrix, i.e., $ \u2211 i = 1 N occ f \phi i f ( r ) \u27e8 \phi i f  \xi ( w ) \u27e9 D f $, in each SCF iteration, while in oefsDFT, the projection of χ⟩ is performed only once at the beginning of the calculation.
VI. SUMMARY
In this work, we have developed an approach to reduce the statistical fluctuations in the electron density, total energy, and forces on the nuclei within the stochastic DFT framework without increasing the number of stochastic orbitals. This achievement was made possible by combining the overlapped embeddedfragmented stochastic DFT^{27} with the energy window stochastic DFT.^{32} The new approach builds on both realspace and energyspace fragmentation, resulting in a significant reduction in the noise in singleparticle observables without affecting the computational time. The performance of the ewefsDFT was tested for a Gcenter embedded in bulk silicon with a small fundamental gap and ingap impurity states, making this a rather challenging system for DFT. Compared to oefsDFT and ewDFT (not shown explicitly here), the statistical error in the forces is ∼6 times smaller in ewefsDFT, resulting in a reduction of ≈30 in the computational wall time. This reduction in noise/computational time is important to accurately describe the structural properties of extended systems without the need to increase the number of stochastic orbitals. Application of the ewefsDFT method to structural minimization is currently under way.^{33}
SUPPLEMENTARY MATERIAL
The asymptotic behavior of singleparticle observables and analysis of the variance of the electron density for ewefsDFT are provided in the supplementary material.
ACKNOWLEDGMENTS
We acknowledge support from the Center for Computational Study of Excited State Phenomena in Energy Materials (C2SEPEM) at the Lawrence Berkeley National Laboratory, which is funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division under Contract No. DEAC0205CH11231 as part of the Computational Materials Sciences Program. Computational resources were provided by the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DEAC0205CH11231. R.B. gratefully acknowledges support from the GermanyIsrael Foundation (GIF) (Grant No. I26303.22018).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.