Linear scaling density functional theory is important for understanding electronic structure properties of nanometer scale systems. Recently developed stochastic density functional theory can achieve linear or even sublinear scaling for various electronic properties without relying on the sparsity of the density matrix. The basic idea relies on projecting stochastic orbitals onto the occupied space by expanding the Fermi-Dirac operator and repeating this for Nχ stochastic orbitals. Often, a large number of stochastic orbitals are required to reduce the statistical fluctuations (which scale as Nχ1/2) below a tolerable threshold. In this work, we introduce a new stochastic density functional theory that can efficiently reduce the statistical fluctuations for certain observable and can also be integrated with an embedded fragmentation scheme. The approach is based on dividing the occupied space into energy windows and projecting the stochastic orbitals with a single expansion onto all windows simultaneously. This allows for a significant reduction of the noise as illustrated for bulk silicon with a large supercell. We also provide theoretical analysis to rationalize why the noise can be reduced only for a certain class of ground state properties, such as the forces and electron density.

Accurate modeling of the electronic structure is important for understanding the thermal, optical, and magnetic properties of molecules and condensed matters. Density functional theory (DFT)1 within the Kohn-Sham (KS) scheme2 is often the method of choice due to its relatively low computational cost, which formally scales as Ne3, where Ne is the number of electrons in the system. However, even this mild power-law scaling limits the applicability of DFT to study complex systems that require modeling of more than 104 electrons.3–5 In order to reduce the scaling, much effort has been devoted to the development of linear scaling DFT algorithms.6–17 The majority of these algorithms either divide a system into subsystems14–17 and attempt to solve the electronic structure problem in each subsystem or rely on the sparsity of the density matrix.6,10,12 Although linear scaling DFT has achieved many successes, difficulties exist in both approaches. Algorithms based on solving subsystem such as like embedding methods require a proper treatment of boundaries between subsystems16,18–22 while methods based on the locality of the density matrix often suffer from slow convergence when the fundamental band gap is small.

Recently, we have introduced an alternative scheme for linear-scaling DFT based on the notion of stochastic orbitals, termed stochastic DFT (sDFT).23 The idea is to compute the electron density and ground state electronic properties with stochastic orbitals that are random linear combinations of the deterministic occupied orbitals. Importantly, the sDFT approach depends on neither the sparsity of density matrix nor how the boundaries between fragments in an embedding scheme are treated. This has been shown for a variety of systems, including insulators and semiconductors.23–28 Although sDFT is very efficient for large systems, lowering the statistical noise of ground state electronic properties to below tolerable thresholds often requires hundreds to thousands of stochastic orbitals, leading to significant computational costs.23 Various techniques based on embedded fragmentation algorithms have been successfully introduced to reduce the noise in the electron density for specific systems, such as molecular systems24 and covalently bonded materials,25,26 but more robust noise reduction schemes are needed, particularly for systems for which it is nontrivial to construct embedded fragments.

Here, we propose a general improvement to sDFT that reduces the noise in the electron density and nuclear forces by breaking the electron density into components corresponding to orbital energy windows and using sDFT to sample each window. We refer to this new method as “energy window stochastic density functional theory” (ew-sDFT). We find that ew-sDFT successfully reduces the statistical noise in the electron density and in certain other ground state electronic structure properties (e.g., forces).

This paper is organized as follows: We briefly review the sDFT approach and then present the energy-window modification followed by a numerical illustration for bulk silicon with a large supercell. Finally, we provide a discussion of the noise reduction mechanism in ew-sDFT.

We begin by briefly reviewing the stochastic density functional theory (for more details, consult Refs. 23–26 and 28). Without loss of generality, we limit the discussion to a system in a box of volume V with periodic boundary condition. Electronic properties are calculated by KS-DFT with Hamiltonian ĥKS given by

ĥKS=t^+v^nl+v^loc+v^h[ρ]+v^xc[ρ],
(1)

where t^, v^nl, v^loc, v^h[ρ], and v^xc[ρ] are kinetic, nonlocal pseudopotential, local pseudopotential, Hartree, and exchange correlation terms. In KS-DFT, the Hartree and exchange correlation potentials explicitly depend on the electron density ρr, which is given by ρ(r)=2i=1Nocc|ψi(r)|2 on Ng real-space grid points, where ψi(r) are the KS orbitals of ĥKS and Nocc is the number of occupied orbitals. ρ(r) can also be expressed as the trace of an operator, i.e.,

ρ(r)=2Tr ρ^δ(rr^),
(2)

where ρ^=limβθβ(ĥKS;μ) is the one-body density matrix, θβ(x; μ) = 1/(1 + eβ(xμ)) is the Fermi-Dirac distribution function with inverse temperature β and chemical potential μ, and δ(rr^) is the Dirac-delta function. In practice, β is chosen to be large enough to converge the ground state properties and μ is tuned such that Tr ρ^=Ne/2 is equal to the number of electron pairs.

In sDFT, ρ(r) in Eq. (2) is determined by using a set of Nχ stochastic orbitals χ (rather than finding all KS orbitals),

ρr=2χρ^δrr^χχ=2ξr2χ,
(3)

where ξ=ρ^χ and ⟨⋯⟩χ imply averaging over an ensemble of stochastic orbitals. One possible choice of χ in real space is χ(r)=±1/ΔV with equal random probabilities of positive or negative values for each grid point, and ΔV = V/Ng represents the volume element. In practice, a set of Nχ stochastic orbitals {χi} is used in sDFT and a sample average leads to ρ¯χ(r)=2Nχi=1Nχ|ξ(r)|2. ρ¯χ(r) is then used to evaluate quantities that are functionals of the density, such as the local pseudopotential energy, Hartree energy, exchange-correlation energy, and the local pseudopotential nuclei force. A ground state quantity corresponding to a one-body operator, Ô, can also be obtained via a trace formula,

O=2Tr(ρ^Ô)=2ξ|Ô|ξχ.
(4)

Equation (4) is used to evaluate kinetic energy (EK) and nonlocal pseudopotential energy (Enl) and nonlocal pseudopotential nuclei force. Applying Eq. (4) with a finite number of samples results in a statistical approximation of the expectation value of Ô.

Achieving linear (or sublinear) scaling in sDFT relies on two requirements. First, applying ρ^=θβ(ĥKS;μ) to a stochastic orbital χ should scale linearly with respect to system size. This is achieved by expanding the Fermi-Dirac operator in a Chebyshev polynomial,29,30

ρ^χ=θβ(ĥKS;β,μ)χ=n=0Ncan(μ,β)Tn(ĥKS)χ,
(5)

where an(μ, β) are the expansion coefficients of the Fermi-Dirac operator31 and Nc is the highest degree of polynomials. Linear scaling then relies on the sparsity of the KS Hamiltonian and not on the density matrix. Second, for many local observables like the energy per particle and the forces on the atoms, the statistical noise does not increase with increasing system size and often decreases with the system size, leading to sublinear scaling.23–26,28

Within the sDFT formulation, the exact electron density ρ(r) can only be recovered by averaging infinitely many stochastic orbitals. Thus, for a finite set of stochastic orbitals, ρ¯χ(r) contains a statistical noise. With the choice of random state χ described above, the variance of density calculated from one stochastic orbital, ρχr=2χθβĥKS;μδrr^χ, is approximately given by25 

Varχ(ρχ(r))=2ρ2(r),
(6)

where Varχ(⋯) represents variance in ρ(r). From the central limit theorem, the variance of the electron density ρ¯χ(r) evaluated with Nχ stochastic orbitals is 2ρ2(r)/Nχ. Equation (6) indicates that the standard deviation of the density is proportional to the density itself. Thus, by dividing ρ(r) into multiple components and using sDFT to sample each of them, it is possible to reduce the statistical noise, as discussed in Sec. III.

The central idea behind the energy window stochastic density functional theory (ew-sDFT) is to divide the space spanned by the occupied orbitals into subspaces and perform sDFT calculation in each subspace before combining contributions from all subspaces. We will illustrate ew-sDFT with two energy windows; however, generalizing the results to multiple windows is straightforward. First, we divide the occupied orbitals, {ψi}, into two subsets: one with orbital energy (εi) smaller than a given energy ε* while orbitals in the other subset have orbital energies higher than ε*, where ε* < μ. The orbital indexes in the two subsets are w1 = {i|εiε*} and w2 = {i|ε* < εi < μ}. The corresponding projection operators onto each energy window are ρ^1=θβ(ĥKS;ϵ*) and ρ^2=θβ(ĥKS;μ)θβ(ĥKS;ϵ*). Similar to sDFT, applying ρ^1 and ρ^2 to a stochastic orbital χ generates filtered stochastic orbitals ξ(1)=ρ^1χ and ξ(2)=ρ^2χ. The deterministic density ρ(r) can be recovered by averaging 2|ξ(1)(r)|2 + 2|ξ(2)(r)|2, i.e.,

ρ(r)=2|ξ(1)(r)|2+|ξ(2)(r)|2χ.
(7)

We evaluate both ρ1^ and ρ2^ using a Chebyshev polynomial expansion given by Eq. (5). We emphasize that ew-sDFT barely increases the computational wall time compared to sDFT because the same Chebyshev polynomials of ĥKS can be used to simultaneously project a random stochastic orbital on both energy windows, generating |ξ(1)⟩ and |ξ(2)⟩ with a similar effort to one window sDFT.

Within ew-sDFT, a general ground state observable can be evaluated according to

O=2Tr(ρ^Ô)=2ξ(1)|Ô|ξ(1)+ξ(2)|Ô|ξ(2)χ
(8)

and by ρχ(r) = 2|ξ(1)(r)|2 + 2|ξ(2)(r)|2, the variance of ρχ(r) in ew-sDFT is given by

Varχ(ρχ(r))=2ρ12(r)+ρ22(r),
(9)

where ρ1r=2iw1|ψir2| and ρ2r=2iw2|ψir2| are densities corresponding to the two energy windows. Recalling that the variance of ρ¯χ(r) from sDFT is 2ρ2(r)=2(ρ1(r)+ρ2(r))2, the variance is thus reduced by 4ρ1(r)ρ2(r). Therefore, with similar computational wall time, ew-sDFT achieves smaller fluctuations in the electron density compared to sDFT.

The implementation of ew-sDFT is similar to sDFT and can be implemented in plane waves, real space grids, or Gaussian basis sets. Starting from the initial guess of the electron density, ew-sDFT first generates Nχ random vectors in real space as described in Sec. II. Next, the chemical potential μ is determined by solving the following equation:

2Nχi=1Nχχi|f(ĥKS;β,μ)|χi=Ne.
(10)

Note that Eq. (10) is the same as the one used in sDFT since χi|ρ^1|χi+χi|ρ^2|χi=χi|ρ^|χi for each |χi⟩. After determining μ, each |χi⟩ is filtered by ρ^1 and ρ^2 using a Chebyshev polynomial expansion as described by Eq. (5) to generate |ξi(1) and |ξi(2). The filtered stochastic orbitals are then used to calculate the density, ρ¯χ(r), and other electronic properties via Eqs. (7) and (8). The KS Hamiltonian is then updated with ρ¯χ(r), and the above calculations repeat until ρ¯χ(r) converges. The same random number seed is used throughout the self-consistency iterations to ensure convergence.

We assessed the performance of ew-sDFT on silicon crystal with various supercell size. Silicon crystal is a challenging system for linear scaling DFT algorithms due to its small band gap within the local density approximation (LDA) to DFT (fundamental gap of ∼0.6 eV). We used β ≈ 600 in inverse Hartree, sufficient to converge the ground state properties with respect to the electron temperature. The same value of β ≈ 600 was also sufficient to converge the ground state properties within ew-sDFT, regardless of the number of energy windows. The wave function cutoff was 30 Ry and the density cutoff was 60 Ry, which corresponds to real space grid spacing of 0.22 Å. Troullier-Martins norm conserving pseudopotentials32 were used with Kleinman-Bylander separable form.33 The application of the nonlocal pseudopotential was performed in real space to reduce the computational costs.34 We first tested a rather small supercell containing 64 silicon atoms as an illustrative example of ew-sDFT. 512 stochastic orbitals were used in both sDFT and ew-sDFT and 10 energy windows were used in ew-sDFT. All statistical quantities were evaluated from 10 independent runs.

In addition to the perfect crystalline silicon structures, we have also compared the performance of sDFT and ew-sDFT on a thermally equilibrated silicon structure. The structure was obtained by equilibrating Si512 crystal at 1000 K and ambient pressure with molecular dynamics simulation using the Stillinger and Weber force-field.35 Due to smaller fundamental gap of thermal equilibrium silicon, we used β ≈ 800 in inverse Hartree to converge the ground state properties in both sDFT and ew-sDFT calculations. All the other computational details were not changed, and all statistical quantities were averaged over 5 independent runs.

In Fig. 1, we plot the standard deviation of ρ¯χ(r) for sDFT (blue curves) and ew-sDFT (red curves). In addition, we plot the estimate of the standard deviation given by Eqs. (6) and (9) for sDFT and ew-sDFT (dashed lines), respectively. We find that the standard deviation of the electron density is roughly proportional to the value of ρ(r), consistent with the theoretical predictions. Both theory and numerical calculations demonstrate that ew-sDFT reduces the statistical noise in the density. With 10 windows, the density fluctuation from ew-sDFT is reduced by nearly a factor of 3 compared to sDFT. Since the standard deviation scales as Nχ1/2, it implies that the sDFT requires about 10 times more stochastic orbitals to achieve the same level of statistical noise in the density.

FIG. 1.

Standard deviation of densities at selected grid points are evaluated from sDFT (solid blue) and ew-sDFT with 10 windows (solid red) calculations for Si64 system. Theoretical estimated density standard deviation from Eqs. (6) and (9) are plotted as dashed blue line for sDFT and dashed red line for ew-sDFT.

FIG. 1.

Standard deviation of densities at selected grid points are evaluated from sDFT (solid blue) and ew-sDFT with 10 windows (solid red) calculations for Si64 system. Theoretical estimated density standard deviation from Eqs. (6) and (9) are plotted as dashed blue line for sDFT and dashed red line for ew-sDFT.

Close modal

A similar picture emerges for the density of states (DOS), as shown in Fig. 2. We constructed the KS Hamiltonian ĥKS[ρ¯χ(r)] from ρ¯χ(r) and diagonalized it in the subspace of occupied orbitals. The DOS of ĥKS[ρ¯χ(r)] withρ¯χ(r) from ew-sDFT is in better agreement with the deterministic result (upper panel of Fig. 2) and the fluctuations in the DOS are smaller (lower panel of Fig. 2).

FIG. 2.

(a) Density of states (DOS) of Si64 by using densities calculated from sDFT (blue line) and ew-sDFT with 10 windows (red line) are compared to deterministic DOS (green line). DOS from sDFT and ew-sDFT are averaged from 10 independent runs. (b) Standard deviations of DOS are shown as blue line for sDFT and red line for ew-sDFT.

FIG. 2.

(a) Density of states (DOS) of Si64 by using densities calculated from sDFT (blue line) and ew-sDFT with 10 windows (red line) are compared to deterministic DOS (green line). DOS from sDFT and ew-sDFT are averaged from 10 independent runs. (b) Standard deviations of DOS are shown as blue line for sDFT and red line for ew-sDFT.

Close modal

Reducing the noise in the electron density should lead to a reduction of the noise of all ground state properties. To test this hypothesis, we calculated the energy per electron and nuclei forces for Si64 and Si512 supercells at T = 0 K and for Si512 supercell at T = 1000 K and compared the results of ew-sDFT to sDFT. 10 energy windows were used for all systems. In addition, to explore the effect of the number of windows, we have used 100 energy windows for the crystal Si512. The energies per electron calculated by sDFT and ew-sDFT agree well with deterministic results, as presented in Table I. Surprisingly, the standard deviations are similar comparing sDFT and ew-sDFT. Increasing the number of windows to 100 does not change this picture.

TABLE I.

The average values of energies per electron, including kinetic energy (Ek), nonlocal pseudopotential energy (Enl), external energy/local pseudopotential energy (Eloc), the summation of Hartree (Eh) and exchange-correlation (Exc) energy and total energy (Etot), are presented in the table. The standard deviations in the last digits of energies are listed in parentheses. Nw is the of energy windows used for ew-sDFT calculation.

SystemMethodNwEkEnlElocEh + ExcEtot
Crystal Si64 sDFT  10.450(40) 6.081(8) −8.633(22) −6.270(6) −26.937(37) 
 ew-sDFT 10 10.450(40) 6.081(8) −8.641(19) −6.271(6) −26.946(36) 
 Deterministic  10.450 6.083 −8.660 −6.266 −26.958 
Crystal Si512 sDFT  10.364(7) 6.105(3) −8.581(6) −6.291(2) −26.969(7) 
 ew-sDFT 10 10.363(6) 6.104(3) −8.587(5) −6.293(2) −26.978(7) 
 ew-sDFT 100 10.363(7) 6.104(3) −8.589(5) −6.293(2) −26.980(7) 
 Deterministic  10.360 6.105 −8.594 −6.292 −26.985 
Thermal Si512 sDFT  10.301(5) 6.104(3) −8.790(5) −6.190(1) −26.942(8) 
 ew-sDFT 10 10.298(6) 6.104(2) −8.800(2) −6.190(1) −26.956(6) 
 Deterministic  10.298 6.104 −8.804 −6.190 −26.959 
SystemMethodNwEkEnlElocEh + ExcEtot
Crystal Si64 sDFT  10.450(40) 6.081(8) −8.633(22) −6.270(6) −26.937(37) 
 ew-sDFT 10 10.450(40) 6.081(8) −8.641(19) −6.271(6) −26.946(36) 
 Deterministic  10.450 6.083 −8.660 −6.266 −26.958 
Crystal Si512 sDFT  10.364(7) 6.105(3) −8.581(6) −6.291(2) −26.969(7) 
 ew-sDFT 10 10.363(6) 6.104(3) −8.587(5) −6.293(2) −26.978(7) 
 ew-sDFT 100 10.363(7) 6.104(3) −8.589(5) −6.293(2) −26.980(7) 
 Deterministic  10.360 6.105 −8.594 −6.292 −26.985 
Thermal Si512 sDFT  10.301(5) 6.104(3) −8.790(5) −6.190(1) −26.942(8) 
 ew-sDFT 10 10.298(6) 6.104(2) −8.800(2) −6.190(1) −26.956(6) 
 Deterministic  10.298 6.104 −8.804 −6.190 −26.959 

Contrary to the case of the total energy, the standard deviation in the nuclei forces is significantly reduced by the introduction of energy windows, as shown in Fig. 3. Increasing the number of windows in ew-sDFT to 100 further reduces the noise in the nuclei forces [Fig. 3 panel (d)], however, the standard deviation seems to converge to a lower limit. This suggests that there should exist an optimal window size (close to 10) balancing between the small additional computational overhead associated with increasing the number of windows and the reduction in the statistical noise. In addition to the role of the number of energy windows, Fig. 3 panel (d) clearly demonstrates that the noise in the nuclei forces does not depend on the system size, implying that ew-sDFT scales linearly with the system size for nuclei forces.

FIG. 3.

Selected average nuclei forces along x axis and standard deviations (error-bars) are shown for (a) crystal Si64, (b) crystal Si512, and (c) thermal Si512. Panel (a): Blue and green symbols are nuclei forces from sDFT and ew-sDFT with 10 windows, respectively. Panel (b): Blue, green, and red dots are nuclei forces from sDFT, ew-sDFT with 10 windows, and ew-sDFT with 100 windows, respectively. Panel (c): The nuclei forces from for a thermally equilibrated Si512 structure using sDFT (blue symbols) and ew-sDFT (green symbols). Since the forces are significantly larger than in the perfectly ordered case, we subtract the corresponding deterministic value from each force. (d) The standard deviations of nuclei forces along the x axis averaged over all nuclei as a function of number of energy windows. The results for Si64 crystal, Si512 crystal, and thermally equilibrated Si512 are shown as blue, red, and green symbols, respectively. sDFT results correspond to 1 energy window.

FIG. 3.

Selected average nuclei forces along x axis and standard deviations (error-bars) are shown for (a) crystal Si64, (b) crystal Si512, and (c) thermal Si512. Panel (a): Blue and green symbols are nuclei forces from sDFT and ew-sDFT with 10 windows, respectively. Panel (b): Blue, green, and red dots are nuclei forces from sDFT, ew-sDFT with 10 windows, and ew-sDFT with 100 windows, respectively. Panel (c): The nuclei forces from for a thermally equilibrated Si512 structure using sDFT (blue symbols) and ew-sDFT (green symbols). Since the forces are significantly larger than in the perfectly ordered case, we subtract the corresponding deterministic value from each force. (d) The standard deviations of nuclei forces along the x axis averaged over all nuclei as a function of number of energy windows. The results for Si64 crystal, Si512 crystal, and thermally equilibrated Si512 are shown as blue, red, and green symbols, respectively. sDFT results correspond to 1 energy window.

Close modal

Turning to discuss the role of thermal fluctuations and the introduction of disorder, comparing the standard deviation of the nuclei forces in Fig. 3 panels (c) and (d), we find that the reduction in the noise in the nuclei forces using 10 energy windows is not limited to a perfectly ordered structure at T = 0 and a similar reductions are observed for the thermally equilibrated structure at T = 1000 K. This signifies that the sort of self-averaging observed here is not a consequence of the perfectly ordered structure.

In this section, we analyze the fluctuations in observables described by the trace over a generic operator Ô and rationalize the behavior depicted in Table I and Fig. 3, where reduction in the fluctuations was observed for the forces but not the total energies per particle. We begin by focusing on a simple case in which the operator Ô does not depend on the electron density, ρ¯χ(r). In the stochastic formulation [cf. Eq. (4)], the expectation value is approximated by Oχ = 2⟨ξ|Ô|ξ⟩. The variance of Oχ can be calculated by the following equation:

Varχ(Oχ)=8i,j=1Noccψi|Ô|ψj2,
(11)

where ψi and ψj are deterministic KS occupied orbitals. We need to generate the KS orbital for the analysis, but of course, ew-sDFT does not require any knowledge of the KS orbitals. The variance given by Eq. (11) is determined by the Frobenius norm of the matrix representation of Ô in the occupied subspace. In the case of ew-sDFT, the variance is given by the sum of variances for each window (for simplicity we take only two windows)

Varχ(Oχ)=8i,jw1ψi|Ô|ψj2+i,jw2ψi|Ô|ψj2.
(12)

Recalling that w1 (w2) is the set of orbital indices for the first (second) energy window, Eq. (12) suggests that the variance of Oχ is propositional to the Frobenius norms of diagonal blocks of ⟨ψi|Ô|ψj⟩ (see Fig. 4). If ⟨ψi|Ô|ψj⟩ is diagonal dominant, the variance calculated from Eqs. (11) and (12) is nearly identical. In this case, ew-sDFT is less efficient at reducing the statistical fluctuations in Oχ. In the limiting case in which ψi are eigenfunctions of Ô, the noise level of Oχ calculated from ew-sDFT is identical to that of sDFT (a single window ew-sDFT). On the hand, if ⟨ψi|Ô|ψj⟩ is off-diagonal dominant, significant noise reduction can be achieved by ew-sDFT, as suggested by comparing Eqs. (11) and (12).

FIG. 4.

The variance of Oχ depends on (a) the Frobenius norm of the whole matrix (shaded zone) in sDFT and (b) the Frobenius norm of the diagonal blocks (shaded zone) in ew-sDFT.

FIG. 4.

The variance of Oχ depends on (a) the Frobenius norm of the whole matrix (shaded zone) in sDFT and (b) the Frobenius norm of the diagonal blocks (shaded zone) in ew-sDFT.

Close modal

To illustrate this, consider the stochastic fluctuations in the kinetic energy, nonlocal pseudopotential energy, and local pseudopotential energy, which are described by the operators t^, v^nl, and v^loc [independent on the electron density, ρ¯χ(r)]. The corresponding matrix elements calculated in the KS eigenstates representation are shown in Fig. 5 panels (a)–(c), respectively. The matrices are clearly diagonal dominant for all three cases [panels (a)–(c)]. Thus, it is not surprising that ew-sDFT does not reduce the noise for the corresponding observables, Ek, Enl, and Eloc (as discussed above, see also Table I). On the contrary, matrices of nuclei forces are not diagonal dominant, as illustrated in Fig. 5 panels (d)–(f), consistent with the ew-sDFT results shown in Fig. 3. We note in passing that for amorphous materials, off-diagonal terms in the matrices of t^, v^nl, and v^loc become important, as shown in the supplementary material. However, these matrices are still diagonal dominant and thus noise reduction in energy terms are limited, as illustrated in the thermal equilibrium Si512 example in Table I. On the other hand, for these amorphous materials, the structure of the corresponding matrices for the force remain nondiagonal dominant and ew-sDFT can achieve excellent noise reduction in atom forces, which is also demonstrated by the thermal equilibrium Si512 example in Fig. 3.

FIG. 5.

The matrices ψiÔψj are plotted for (a) Ô=t^, (b) Ô=v^nl, and (c) Ô=v^loc, (d) Ô=v^locX1, (e) Ô=v^nlX1, and (f) Ô=v^locX1+v^nlX1, where X1 is the x coordinate of the first atom.

FIG. 5.

The matrices ψiÔψj are plotted for (a) Ô=t^, (b) Ô=v^nl, and (c) Ô=v^loc, (d) Ô=v^locX1, (e) Ô=v^nlX1, and (f) Ô=v^locX1+v^nlX1, where X1 is the x coordinate of the first atom.

Close modal

Next, we consider the case where the operator Ô depends on the electron density, ρ¯χ(r). This is important for observables like the Hartree energy and the exchange-correlation energy required to fully understand the behavior with respect to the number of energy windows observed for the total energy. Consider the variance of a density functional f[ρ¯χ(r)], where ρ¯χ(r) is calculated from Nχ stochastic orbitals in sDFT or ew-sDFT. Approximating ρ¯χ(r) as ρ¯χ(r)ρ(r)(1+α(r)), where α(r) is a small random perturbation of order ϵ where ϵ = 1/Nχ, and retaining terms to order ϵ leads to

f[ρ¯χ(r)]=f[ρr]+ϵδ2fδρ(r)δρ(r)ρ2(r,r)drdr+O(ϵ3/2),
(13)
Var(f[ρ¯χ(r)])=8ϵi,j=1Noccψi|δfδρ|ψj2+O(ϵ3/2),
(14)

where, as before, ρ(r,r)=2r|ρ^|r is the single particle density matrix. For linear functional like the local pseudopotential energy, f[ρ¯χ(r)]=f[ρr] and the functional variance given by Eq. (14) is identical to applying the central limit theorem to the variance given by Eq. (11) with Ô replaced by a local effective operator δfδρ[ρr](r). For a nonlinear functional, δ2fδρ(r)δρ(r)ρ2(r,r)drdr is nonzero, corresponding to “bias” which decays as 1/Nχ as recently discussed by Cytter et al.27 and Fabian et al.28 In ew-sDFT, the average and variance of f[ρ¯χ(r)] becomes

f[ρ¯χ(r)]=f[ρ(r)]+ϵδ2fδρ(r)δρ(r)ρ12(r,r)+ρ22(r,r)drdr+O(ϵ3/2),
(15)
Var(f[ρ¯χ(r)])=8ϵi,jw1ψi|δfδρ|ψj2+i,jw2ψi|δfδρ|ψj2+O(ϵ3/2).
(16)

Using the above relations in Eqs. (11)–(16) results in the following expression for the variance of total energy:

Var(Etot)=8ϵi,j=1Noccψi|ĥKS|ψj2=8ϵi=1Noccεi2+O(ϵ3/2)
(17)

for sDFT and

Var(Etot)=8ϵi,jw1ψi|ĥKS|ψj2+i,jw2ψi|ĥKS|ψj2+O(ϵ3/2)=8ϵi=1Noccεi2+O(ϵ3/2)
(18)

for ew-sDFT. The expressions for the total energy variance given by Eqs. (17) and (18) are identical, suggesting that the statistical noise in ew-sDFT is identical to that of sDFT, consistent with the results shown above. We wish to emphasize that the results in Eqs. (17) and (18) also support the observed sublinear scaling of the energy per electron with both sDFT and ew-sDFT, since the occupied orbital energies are bounded, and thus, the variance of the total energy is proportional to the system size. This implies that the standard deviation of the energy per electron is inverse proportional to the square root of the system size, explaining our previous observations.23,24

In this study, we developed a stochastic density functional theory scheme which allows us to significantly reduce the statistical noise in certain ground state observable with DFT. This was achieved by dividing the occupied space into energy windows and using correlated sampling in each window. Our new scheme (ew-sDFT) is able to reduce the statistical noise in both the electron density and nuclei forces without significantly increasing the computational time, but not for the total energy or total energy per electron. We showed that noise reduction of an observable in ew-sDFT correlates with whether or not the corresponding operator is diagonal dominant in KS eigenstate representation. The total energy is diagonal dominant, and no noise reduction was observed. However, for the forces and electron density, this was not the case and we observed a significant reduction in the noise that results in nearly an order of magnitude reduction in the computational effort compared to sDFT (since the noise is proportional to 1/Nχ). This reduction allows for a more rapid sampling of the canonical distribution using methods that do not rely on the total energy, such as Langevin dynamics.26 The developed ew-sDFT provides an additional approach to noise reduction within the family of sDFT methods and can be combined with embedded fragmentation schemes to further improve the accuracy and efficiency of stochastic density functional theories.

Detailed derivations for variance of ground state properties in sDFT and ew-sDFT are provided in the supplementary material. The supplementary material also includes operator matrices for disordered Si64 system.

R.B. gratefully acknowledges support from the Israel Science Foundation, Grant No. 800/19. D.N. and E.R. are grateful for support by 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. DEAC02-05CH11231 as part of the Computational Materials Sciences Program. Resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231 are greatly acknowledged.

1.
P.
Hohenberg
and
W.
Kohn
,
Phys. Rev.
136
,
B864
(
1964
).
2.
W.
Kohn
and
L. J.
Sham
,
Phys. Rev.
140
,
A1133
(
1965
).
3.
S.
Baroni
,
S.
de Gironcoli
,
A.
Dal Corso
, and
P.
Giannozzi
,
Rev. Mod. Phys.
73
,
515
(
2001
).
4.
T.
Frauenheim
,
G.
Seifert
,
M.
Elstner
,
T.
Niehaus
,
C.
Kohler
,
M.
Amkreutz
,
M.
Sternberg
,
Z.
Hajnal
,
A. D.
Carlo
, and
S.
Suhai
,
J. Phys.: Condens. Matter
14
,
3015
(
2002
).
5.
C.
Freysoldt
,
B.
Grabowski
,
T.
Hickel
,
J.
Neugebauer
,
G.
Kresse
,
A.
Janotti
, and
C. G.
Van de Walle
,
Rev. Mod. Phys.
86
,
253
(
2014
).
6.
F.
Mauri
,
G.
Galli
, and
R.
Car
,
Phys. Rev. B
47
,
9973
(
1993
).
7.
P.
Ordejón
,
D. A.
Drabold
,
M. P.
Grumbach
, and
R. M.
Martin
,
Phys. Rev. B
48
,
14646
(
1993
).
8.
S.
Goedecker
,
J. Comput. Phys.
118
,
261
(
1995
).
9.
E.
Hernández
and
M. J.
Gillan
,
Phys. Rev. B
51
,
10157
(
1995
).
11.
A. H. R.
Palser
and
D. E.
Manolopoulos
,
Phys. Rev. B
58
,
12704
(
1998
).
12.
R.
Baer
and
M.
Head-Gordon
,
Phys. Rev. Lett.
79
,
3962
(
1997
).
13.
S.
Goedecker
,
Rev. Mod. Phys.
71
,
1085
(
1999
).
15.
16.
T.
Zhu
,
W.
Pan
, and
W.
Yang
,
Phys. Rev. B
53
,
12713
(
1996
).
17.
F.
Shimojo
,
R. K.
Kalia
,
A.
Nakano
, and
P.
Vashishta
,
Phys. Rev. B
77
,
085103
(
2008
).
18.
W.
Yang
and
T.
Lee
,
J. Chem. Phys.
103
,
5674
(
1995
).
19.
A. W.
Götz
,
S. M.
Beyhan
, and
L.
Visscher
,
J. Chem. Theory Comput.
5
,
3161
(
2009
).
20.
T. A.
Wesolowski
,
S.
Shedge
, and
X.
Zhou
,
Chem. Rev.
115
,
5891
(
2015
).
21.
J. D.
Goodpaster
,
T. A.
Barnes
, and
T. F.
Miller
,
J. Chem. Phys.
134
,
164108
(
2011
).
22.
C.
Huang
and
E. A.
Carter
,
J. Chem. Phys.
135
,
194104
(
2011
).
23.
R.
Baer
,
D.
Neuhauser
, and
E.
Rabani
,
Phys. Rev. Lett.
111
,
106402
(
2013
).
24.
D.
Neuhauser
,
R.
Baer
, and
E.
Rabani
,
J. Chem. Phys.
141
,
041102
(
2014
).
25.
M.
Chen
,
R.
Baer
,
D.
Neuhauser
, and
E.
Rabani
,
J. Chem. Phys.
150
,
034106
(
2019
).
26.
E.
Arnon
,
E.
Rabani
,
D.
Neuhauser
, and
R.
Baer
,
J. Chem. Phys.
146
,
224111
(
2017
).
27.
Y.
Cytter
,
E.
Rabani
,
D.
Neuhauser
, and
R.
Baer
,
Phys. Rev. B
97
,
115207
(
2018
).
28.
M. D.
Fabian
,
B.
Shpiro
,
E.
Rabani
,
D.
Neuhauser
, and
R.
Baer
, “
Stochastic density functional theory
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
(published online).
29.
R.
Kosloff
,
J. Phys. Chem.
92
,
2087
(
1988
).
30.
R.
Kosloff
,
Annu. Rev. Phys. Chem.
45
,
145
(
1994
).
31.
R.
Baer
and
M.
Head-Gordon
,
J. Chem. Phys.
107
,
10003
(
1997
).
32.
N.
Troullier
and
J. L.
Martins
,
Phys. Rev. B
43
,
1993
(
1991
).
33.
L.
Kleinman
and
D. M.
Bylander
,
Phys. Rev. Lett.
48
,
1425
(
1982
).
34.
R. D.
King-Smith
,
M. C.
Payne
, and
J. S.
Lin
,
Phys. Rev. B
44
,
13063
(
1991
).
35.
F. H.
Stillinger
and
T. A.
Weber
,
Phys. Rev. B
31
,
5262
(
1985
).

Supplementary Material