{'fraction_non_alphanumeric': 0.028930817610062894, 'fraction_numerical': 0.006415094339622641, 'mean_word_length': 4.7280259365994235, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Self-averaging stochastic Kohn-Sham density functional theory

We formulate the Kohn-Sham density functional theory (KS-DFT) as a statistical theory in which the electron density is deter-mined from an average of correlated stochastic densities in a trace formula. The key idea is that it is sufficient to converge the total energy per electron to within a predefined statistical error in order to obtain reliable estimates of the electronic band structure, the forces on nuclei, the density and its moments, etc. The fluctuations in the total energy per electron are guaranteed to decay to zero as the system size increases. This facilitates"self-averaging"which leads to the first ever report of sublinear scaling KS-DFT electronic structure. The approach sidesteps calculation of the density matrix and thus is insensitive to its evasive sparseness, as demonstrated here for silicon nanocrystals. The formalism is not only appealing in terms of its promise to far push the limits of application of KS-DFT, but also represents a cognitive change in the way we think of electronic structure calculations as this stochastic theory seamlessly converges to the thermodynamic limit.

A complimentary DFT formulation is based on orbital-free methods [17][18][19], which promise linear-scaling in algorithmic complexity regardless of the dimensionality of the system. However, the imperfect representation of the kinetic energy in these approaches may hamper their accuracy. A recent promising development along these lines, based on the availability of an approximate functional connecting the ground state density of non-interacting electrons to the underlying potential, provides an estimate of the kinetic energy from adiabatic connection integration [20].
In this letter we develop a theory that allows an in-principle accurate bridge between the Kohn-Sham and the orbital-free approaches to DFT, gaining from both worlds. The core idea is to use a stochastic technique to DFT (SDFT) calculating directly the density from the KS Hamiltonian using a trace formula without computing KS orbitals or density matrices. The use of a stochastic technique combined with DFT has been put forward several decades ago based on path integral Monte Carlo [21][22][23][24], but has not found general use due to algorithmic problems [21]. In contrast, our approach uses a deterministic Chebyshev expansion of the density matrix projection operator and applies it to stochastic orbitals thereby circumventing the pathologies of path integral Monte Carlo.
The stochastic aspect guarantees an unbiased statistical error distribution and thus enables us to demand convergence, not in the total energy as typically required, but in the total energy per electron. Statistical mechanics guarantees that in the thermodynamic limit the fluctuations in the total energy per electron are negligible. This expected decrease of fluctuations is the basis for the new concept of "self-averaging" introduced in the present electronic structure theory. It also enables the development of an algorithm with sub-linear scaling in computational complexity. Thus, the present approach provides for a smooth transition from traditional quantum methods originally developed for small systems to the realm of macroscopic structures understood in terms of ensembles and statistical fluctuations.
Related approaches have been recently developed by us for estimating the rate of multiexciton generation in nanocrystals (NCs) [25], for a linear scaling calculation of the exchange energy [26], for overcoming the computational bottleneck in Møller-Plesset second order perturbation theory [27], and for calculating the Random-Phase-Approximation correlation energy for DFT [28]. These are corrections to the fundamental problem of describing the electron density and total energy, which is addressed here for the first time.
In DFT, given a potential ( ), the electron density can be calculated as a trace: where ( ) [ ] is a smoothed Heaviside step function with , is the frontier orbital gap, ̂ ( ̂) is the one-body Hamiltonian (atomic units are used throughout), and ̂( ) ( ̂). The number of electrons is determined by the chemical potential parameter and its value ( ) is chosen to yield the desired number of electron in the system ( ): In the usual approach to KS-DFT the trace is performed using the lowest eigenstates ( ) of ̂, where ( ̂) projects on the occupied space so the density is given by ( ) ∑ | ( )| . The cubic scaling of conventional KS-DFT methods results from the need to find the lowest energy eigenstates of the Hamiltonian.
In order to eliminate the need for determining the KS eigenstates we use stochastic wave-functions, ( ), for which the completeness property holds: where ̂ is the identity operator and ⟨ ⟩ represents a statistical expectation value (mean) over the stochastic wavefunctions. From equations (1) and (3) the density is given by: where ( ) is a stochastic occupied orbital evaluated using a Chebyshev expansion [15,[29][30][31]: In the above, is the length of the Chebyshev polynomial, | ⟩ is a fixed linear combination of | ⟩, ̂| ⟩ and | ⟩, and (see details in Ref. [32]). The Chebyshev expansion length is proportional to , where is the energy range of the Hamiltonian ̂ and is of the order . Note that while most of the numerical effort goes into the calculation of , these states do not depend on the chemical potential . Only the coefficients ( ) depend on it. Thus it is not significantly more expensive to determine | ⟩ for several values of than it is for just one. Stochastic wave functions have been used before for estimating the density of states [33] but the present work is the first to use them for computing the electronic spatial density with an SCF KS-DFT procedure.
Since ⟨ ( ) ( ) ⟩ is the KS density matrix, any oneparticle observable (e.g. the kinetic energy, the local and nonlocal pseudopotential energy, forces on the nuclei, etc.) is computed as an average: In the SDFT calculation for electrons, one starts by guessing an initial density ( ) and the chemical potential . For example, the initial density can be taken as the sum over the spherical densities of the atoms. The following steps are then used to update the density and the chemical potential:  Compute the KS potential from the density.  Generate stochastic orbitals. In a grid (or plane-wave) representation ( ) √ where is a random phase and is the grid spacing. For each orbital, calculate its projection onto the occupied space ( ) using Eq. (5).
 Determine the density ( ) from Eq. (4) for several values of using the same Chebyshev expansion. Linearly interpolate to find the correct value of the chemical potential, using Eq. (2).  Reiterate until converges below a defined threshold.
It is important to note that one must keep the same random orbitals throughout the self-consistent procedure. Otherwise, the convergence of the self-consistent procedure would be limited by the stochastic noise. In addition, while the above procedure is exact in the limit of an infinite set of , in practice, one uses a finite set containing stochastic orbitals and the estimates of the density and the expectation values involve a statistical error (SE) that is proportional to √ . Table 1: Parameters for the silicon nanocrystals. Shown, are the number of grid points in each Cartesian dimension , the parameter used to represent the Heaviside function, the length of the Chebyshev expansion ( ), the total energy per electron from the deterministic LDA calculation (when available) and the corresponding stochastic result (based on iterations).
We demonstrate the approach on a series of silicon nanocrystals (NCs) described in Table 1. The calculations were done using the local density approximation (LDA) for the XC energy [34] and the nuclear electron interaction was described using Troullier-Martins norm conserving pseudopotentials [35] with a local and nonlocal part: ̂ ( ) ̂ . The total energy was converged self-consistently to per electron using DIIS acceleration [36]. The last two columns in Table 1 show several meV deviations of the stochastic energy per particle from the exact LDA value. The expectation value was estimated from 15 samples of independent SCF runs, each with a different set of stochastic orbitals. The SE is the square root of the variance of this random variable and is estimated from the standard deviation of the sample.
In Figure 1 we show the SE as a function of system size. As the system size increases the SE decreases. This can be explained by considering a large system composed of two noninteracting identical parts having and electrons. The SE in the total energy of the two independent parts relates to the SE in the total energy of each part as . Thus, for large systems the SE of the total energy scales as where is the total number of electrons. For the statistical error of the energy per particle one then finds that ( ) , i.e. inversely proportional to the system size. Thus the SE per electron is expected to decrease as for large system. Figure 1 shows that this expected behavior holds for systems with .  Figure 2 shows the CPU time for the SDFT as a function of system size for a fixed SE per electron of 10 meV. We find that the algorithmic scaling is sub-linear with system size. The total CPU time depends on several factors: The grid size which scales linear with the number of electrons; the inverse HOMO-LUMO gap which determines the values of and the length of the Chebyshev expansion and roughly scales as ( is the bulk band gap and is a constant); and the number of stochastic orbitals needed to reach a given SE per electron.
The fact that the SE per electron decreases as √ implies that, for a fixed SE per electron, fewer stochastic orbitals are required to converge the expectation values of the energy per electron to within a desired tolerance as the system grows (see inset of Figure 2). In the limit of a macroscopic system only a single stochastic orbital will be needed to estimate the total energy per electron with negligible error. This "selfaveraging" property leads to sub-linear scaling at intermediate system sizes. Once the system size is so large as to enable a single stochastic orbital to be used (for the present system, this is extrapolated to be at 500,000 electrons), the numerical effort of the approach will scale linearly with system size (regime not shown).
In comparison to the deterministic approach, the CPU time for the SDFT calculation crosses that of the deterministic calculation at systems containing about 3000 electrons (order of 700 silicon atoms) for a SE per electron of 10meV. It is interesting to note that the calculations shown in Figure 1 and Figure 2 were carried on a parallel cluster (with 180 cores) for the stochastic approach and on a single core for the deterministic one. As the communication overhead is negligibly small, such a comparison can be made without any bias towards the single-core calculation. This is another significant advantage of the present theoryit scales linearly with the number of cores as long as this number is smaller than The scaling of the present approach with CPU time and memory is unprecedented for 3D structures. Once again, it relies on converging the total energy per electron rather than the traditional approaches to converge the total energy. The convergence of the total energy per electron is perhaps sufficient to obtain an accurate estimate of this property alone. However, the SE in the total energy itself is rather high, far from a desired accuracy of tens of meVs. This raises a concern whether the current approach can be used as a reliable tool to obtain expectation values of quantities that depend on all electrons. To address this concern, in Figure 3 we plot the forces along an arbitrary ( ) direction obtained by the SDFT for one of the silicon NCs studied. Calculations for the other NCs provide a similar picture. The forces were obtained by integrating the density with the gradient of the electronnucleolus potential.
The two panels in Figure 3 correspond to two sets of SDFT calculations in which a small ( ) and a larger ( ) number of stochastic orbitals were used. The SE in the total energy is ~6eV and ~3eV in the two calculations, respectively, yet the error in the forces is of statistical nature and is perfectly controllable: it drops by a factor 2 when I increases by a factor of 4. The SE in the force on the silicon atoms (low index in the figure) is about is 4 times larger than that on the H atoms (high index in the figure). This is due to the higher number of electrons near the silicon nucleus (also by a factor of 4) and could likely be reduced with importance sampling techniques, as will be explored in future studies. The error bar on the force on a silicon atom when the SE in the total energy per electron is 10meV is ~0.02 such a value is low enough to enable Langevin dynamics simulations which shed light on the equilibrium and non-equilibrium statistical mechanical properties of large systems. Despite the orbital-less nature of the calculation, the density of states (DOS) and KS orbital energies are readily available.
In the left panel of Figure 4 we show the DOS for for energies close to the chemical potential estimated from a finite difference formula applied to ( ) . Results obtained from two different stochastic runs are compared to those of the direct approach using the KS orbital gies ( ) ∑ ( ( )) . The value of the chemical potential, for each stochastic run varies ( typically "sticks" to near to the HOMO or to the LUMO of the system) but the DOS is similar between the different runs and agrees well with that of the direct approach, to within a small statistical error.
In the right panel of Figure 4 we compare the value of the KS orbital energies for the same NC. The agreement between the stochastic and deterministic approach is remarkable. The inset of Figure 4 shows the KS orbital energy error defined by the difference between the stochastic and deterministic estimates near the top of the valance band and the bottom of the conduction band. The small error ( ) has a stochastic and deterministic components. The former can be decreased by increasing the number of stochastic orbitals used for the SDFT procedure while the latter is controlled by the length of the Chebyshev polynomial.
In summary, we have presented a stochastic formulation of DFT that provides a link between KS and orbital-free DFT formulations. The electron density is given in terms of a trace formula which is evaluated using stochastic occupied orbitals generated by a Chebyshev expansion of the occupation operator, rather than by finding all the KS orbitals. Due to the statistical nature of the formalism, the electron density and ex-pectation values involve an unbiased statistical error, which obeys normal large-number statistics and thus, can be controlled by increasing the number of stochastic orbitals . We demonstrated the method and its properties on silicon NCs.
A novel aspect of our formalism is that it allows to control the total energy per electron instead of the total energy itself while obtaining statistically significant estimates of electron densities forces on nuclei and band structure from the calculation. For this reason, the method seamlessly connects with the thermodynamic limit where the statistical fluctuations in the total energy per electron vanish. Therefore, our formalism enjoys an effect of "self-averaging", whereas the silicon NC grows one requires fewer stochastic iterations for reaching the same accuracy. It is this effect which enables sublinear scaling, never before seen in KS-DFT electronic structure.
Finally, since the method produces a KS Hamiltonian, it can be integrated with the MP2 and RPA stochastic approaches we recently developed [27,28].