Stochastic density functional theory (sDFT) is becoming a valuable tool for studying ground-state 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 sub-linear scaling of certain ground-state 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 χ 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 ground-state properties, such as the electron density, total energy, and forces on the nuclei, as demonstrated for a G-center in bulk silicon.

Kohn–Sham (KS) density functional theory1,2 (DFT) is widely used to study a wide range of systems due to its capability of quantitatively predicting ground-state properties at a moderate computational cost of O ( N e 3 ) , where Ne 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 104 electrons or more, such as nanostructures,3 complex materials,4 and large biomolecules,5 is still a severe challenge for today’s DFT implementations. Linear-scaling methods for DFT based on dividing the entire system into subsystems6–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 linear-scaling 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 ground-state properties. In practice, the required number of stochastic orbitals does not increase with the system size for evaluating many ground-state properties,26,27 leading to linear scaling or even sub-linear 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 conditions28–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 embedded-fragmented stochastic density functional theory” (o-efsDFT).27] Recently, we introduced an alternative technique to mitigate the statistical noise, referred to as “energy-window sDFT” (ew-sDFT),32 where the occupied space was divided into energy-resolved 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 embedded-fragmented 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 embedded-fragmented stochastic DFT (ew-efsDFT) 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 ew-sDFT or o-efsDFT, as illustrated for a G-center 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 o-efsDFT and ew-sDFT methods, both central to the development of the current noise reduction scheme. In Sec. IV, we provide the details of the proposed ew-efsDFT and a summary of the algorithm. Assessment of the new approach for a challenging G-center embedded in bulk silicon is presented in Sec. V alongside a discussion of the computational complexity and cost of the ew-efsDFT. Finally, in Sec. VI, we summarize the main developments.

Consider an extended system described by KS-DFT, with a KS Hamiltonian ( h ̂ K S ) given by

h ̂ K S = t ̂ + v ̂ nl + v ̂ loc + v ̂ H [ ρ ] + v ̂ xc [ ρ ] ,
(1)

where t ̂ , v ̂ nl , v ̂ loc , v ̂ H [ ρ ] , and v ̂ xc [ ρ ] are the operators of the kinetic energy, the non-local 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, ρ r , which is formally given by (we assume closed-shell and ignore spin–orbit couplings for simplicity)

ρ ( r ) = 2 Tr ρ ̂ δ ( r r ̂ ) = lim β 2 Tr θ β ( h ̂ KS , μ ) δ ( r r ̂ ) ,
(2)

where ρ ̂ = θ β ( h ̂ KS , μ ) is the one-body 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, Ne = ∫drρ(r). Other smooth functions to approximate a step function can also be used instead of θβ(xμ). In KS-DFT, the electron density can also be written in terms of the KS orbitals (eigenstates of the KS Hamiltonian), ϕi(r),

ρ ( r ) = 2 i N occ | ϕ i ( r ) | 2 ,
(3)

where Nocc is the number of occupied orbitals.

In sDFT, the trace in Eq. (2) is replaced by averaging the expectation value of θ β ( h ̂ KS , μ ) δ ( r r ̂ ) ,

ρ ( r ) = 2 χ θ β ( h ̂ KS , μ ) δ ( r r ̂ ) χ χ ,
(4)

where |χ⟩ is a stochastic orbital and ⟨⋯⟩χ implies averaging over an ensemble of stochastic orbitals. The stochastic orbitals are represented on a real-space grid with Ng grid points; each grid point is assigned a random value ± 1 / Δ V , where ΔV = V/Ng is the volume element and V is the volume of the supercell. Equation (4) can be rewritten in a compact form as

ρ ( r ) = 2 | ξ ( r ) | 2 χ ,
(5)

where |ξ⟩ is a projected stochastic orbital, which is given as

| ξ = ρ ̂ | χ = θ β ( h ̂ KS , μ ) | χ .
(6)

The projection of the stochastic orbitals onto the occupied space is obtained by expanding θ β ( h ̂ KS , μ ) in Chebyshev polynomials,34,35

θ β ( h ̂ KS , μ ) = n = 0 N c a n ( μ , β ) T n ( h ̂ KS ) ,
(7)

where Nc is the length of the Chebyshev polynomial expansion, an(μ, β) are the expansion coefficients, and Tn are the Chebyshev polynomials of order n.

Ground-state observables corresponding to any one-body operator, O ̂ , can be evaluated using a similar stochastic trace formula,

O = 2 Tr ( ρ ̂ O ̂ ) = 2 ξ | O ̂ | ξ χ .
(8)

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 χ 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 103 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.

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 KS-DFT. 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 embedded-fragmented stochastic DFT (o-efsDFT), which is central to the proposed ew-efsDFT. Full details of the approach can be found elsewhere.27 

In o-efsDFT, the supercell is divided into fragments referred to as “core regions” (see Fig. 1 for an illustration) Cf wrapped by “buffer regions” (Bf) to form dressed fragments (Df = CfBf), where f is the fragment index. The fragment density matrix, ρ ̂ f , is given by

r | ρ ̂ f | r = i = 1 N occ f r | φ i f φ i f | r , r D f 0 , r D f
(9)

for rCf. In the above equation, φf(r) are the KS orbitals for fragment f obtained from a deterministic KS-DFT 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:

ρ ( r ) = 2 f r | ρ ̂ f ρ ̂ f | r + 2 | ξ ( r ) | 2 χ 2 f r | ρ ̂ f | χ χ | ρ ̂ f | r χ = 2 f ρ f ( r ) + 2 | ξ ( r ) | 2 χ 2 f | ξ f ( r ) | 2 χ ,
(10)

where the fragment electron density is ρ f ( r ) = i = 1 N occ f | φ i f ( r ) | 2 , ξ f ( r ) = i = 1 N occ f φ i f ( r ) φ i f | χ D f , and φ i f | χ D f = D f d r φ i f ( r ) * χ ( r ) . We use the relationship ρ f ( r ) = r | ρ ̂ f ρ ̂ f | r in Eq. (10) since KS-DFT methods in the 0 K limit are adopted to calculate fragment KS orbitals.

FIG. 1.

An illustration of the overlapped fragmented scheme. Each small solid blue/red square represents a core region Cf. Focusing on the solid red square core region of a fragment, the region within the dashed red square corresponds to the dressed fragment, Df. The region between the solid red square and dashed red square is the buffer region.

FIG. 1.

An illustration of the overlapped fragmented scheme. Each small solid blue/red square represents a core region Cf. Focusing on the solid red square core region of a fragment, the region within the dashed red square corresponds to the dressed fragment, Df. The region between the solid red square and dashed red square is the buffer region.

Close modal

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 fragmentation-based DFT approach.8,36 However, unlike divide-and-conquer based methods, the stochastic approach does not introduce systematic errors resulting from the fragmentation in the limit Nχ.

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 ̂ 1 , , P ̂ N w . Here, the projector P ̂ w is defined as

P ̂ w = θ β ( h ̂ KS , ε w ) θ β ( h ̂ KS , ε w 1 )
(11)

for 1 ≤ wNw, where = ε 0 < ε 1 < < ε N w = μ . In ew-sDFT, the electron density is given by the sum of all projected densities,

ρ ( r ) = 2 w = 1 N w ξ ( w ) ( r ) 2 χ 2 w = 1 N w ρ ( w ) ( r ) ,
(12)

where | ξ ( w ) = P ̂ w | χ is a projected stochastic orbital for window w. |ξ(w)⟩ are calculated simultaneously with a single Chebyshev expansion, i.e.,

ξ ( w ) = P ̂ w | χ = n = 0 N c b n ( w ) ( ε w , ε w 1 ) T n ( h ̂ KS ) | χ ,
(13)

where Nc is the length of the Chebyshev polynomial (chosen sufficiently large to converge the expansion), b n ( w ) ( ε w , ε w 1 ) is the expansion coefficient (depends on the window), and Tn(x) is the Chebyshev polynomial of order n. The variance of the electron density in this scheme is 8 w = 1 N w ρ ( w ) ( r ) 2 , which is smaller than the variance in sDFT given by 8 w = 1 N w ρ ( w ) ( r ) 2 .32 

Combining the energy window approach with the fragmentation approach results in the following expression for the electron density at a grid point r:

ρ ( r ) = 2 f ρ f ( r ) + w = 1 N w 2 r | ρ ̂ P ̂ w | χ χ | P ̂ w ρ ̂ | r χ 2 f r | ρ ̂ f P ̂ w | χ χ | P ̂ w ρ ̂ f | r χ = 2 f ρ f ( r ) + 2 w = 1 N w | ζ ( w ) ( r ) | 2 χ 2 f w = 1 N w | ξ f ( w ) ( r ) | 2 χ .
(14)

In the above equation, the projection operators on the energy windows are defined as

P ̂ w = θ β ( h ̂ KS , ε w ) θ β ( h ̂ KS , ε w 1 ) 1 w < N w P ̂ N w = I ̂ w = 1 N w 1 P ̂ w ,
(15)

where { ε } ε 0 ε N w 1 (ɛ0 = −) define the boundaries of the energy windows. We would like to signify the differences between the projection operators used for ew-sDFT [cf. Eq. (11)] and those used in the current ew-efsDFT approach [cf. Eq. (15)]. In ew-sDFT, w = 1 N w P ̂ w equals the density matrix ρ ̂ , while in ew-efsDFT, it equals the unit operator, I ̂ . Furthermore, in ew-sDFT, the highest energy window is set to ε N w = μ , while in ew-efsDFT, the energy windows are held fixed (but not necessarily equally spaced) for the entire self-consistent procedure and are chosen to be independent of the chemical potential, μ. The use of a fixed set of energy windows simplifies the on-the-fly calculations of chemical potential [see Eq. (18)].

In Eq. (14), the action of ρ ̂ P ̂ w and P ̂ w on |χ⟩ is obtained using a proper Chebyshev series,

| ζ ( w ) = ρ ̂ P ̂ w | χ = n = 0 N c a n ( w ) ( μ , ε w , ε w 1 ) T n ( h ̂ KS ) | χ , | ξ ( w ) = P ̂ w | χ = n = 0 N c b n ( w ) ( ε w , ε w 1 ) T n ( h ̂ KS ) | χ .
(16)

Finally, as before, the density of each fragment is given by ρ f ( r ) = i = 1 N occ f | φ i f ( r ) | 2 and the stochastic projected orbitals for each fragment are given by

ξ f ( w ) ( r ) = i = 1 N occ f φ i f ( r ) φ i f | ξ ( w ) D f .
(17)

For a fixed set of energy windows, ɛw, the chemical potential (μ) can be obtained by solving

N ( μ ) = 2 f C f d r ρ f ( r ) + 2 χ | ρ ̂ ( μ ) | χ χ 2 f w = 1 N w C f d r | ξ f ( w ) ( r ) | 2 χ ,
(18)

where C f d r imply that the integrals are performed in real space in region rCf. In the above equation, χ | ρ ̂ ( μ ) | χ is evaluated by expanding ρ ̂ in a Chebyshev series χ | ρ ̂ ( μ ) | χ = n = 0 N c c n ( μ ) χ | T n ( h ̂ KS ) | χ and the chemical potential is determined by solving for N(μ*) = Ne, where Ne is the total number of electrons in the system.

Similar to the electron density, other ground-state observables such as the kinetic energy,

E k = 2 f i = 1 N occ f φ i f | t ̂ | φ i f C f + 2 w = 1 N w ζ ( w ) | t ̂ | ζ ( w ) χ 2 i = 1 N w f ξ f ( w ) | t ̂ | ξ f ( w ) C f χ ,
(19)

or the non-local pseudopotential energy,

E nl = 2 f I , R I C f i = 1 N occ f φ i f | v ̂ nl I | φ i f D f + 2 w = 1 N w I ζ ( w ) | v ̂ nl I | ζ ( w ) χ 2 w = 1 N w f I , R I C f ξ f ( w ) | v ̂ nl I | ξ f ( w ) C f χ ,
(20)

or the non-local pseudopotential contribution to the forces on the nuclei,

F nl I = 2 i = 1 N occ f φ i f v ̂ nl I R I φ i f D f + 2 w = 1 N w ζ ( w ) v ̂ nl I R I ζ ( w ) χ 2 w = 1 N w ξ f ( w ) v ̂ nl I R I ξ f ( w ) D f χ ,
(21)

is expressed in ew-ofsDFT 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 ew-efsDFT method to reduce the noise in the density, energy, and forces on the nuclei can be summarized as follows:

  1. Generate the KS orbitals { φ i f ( r ) } for each dressed fragment and ρ f ( r ) = i = 1 N o c c f | φ i f ( r ) | 2 using a deterministic DFT. ρ(r) = 2fρf(r) is used as the initial electron density guess.

  2. For each stochastic orbital χ(r) (defined above), calculate and store on the grid the projected stochastic orbital ζ ( w ) ( r ) = r | P ̂ w | χ and also store the Chebyshev moments, χ | T n ( h ̂ KS ) | χ .

  3. For each window and for each stochastic orbital, generate and store on the grid ξ f ( w ) ( r ) = i = 1 N occ f φ i f ( r ) φ i f | ξ ( w ) D f .

  4. Solve for μ* (N(μ*) = Ne) with the regula falsi method using the Chebyshev moments ( χ | T n ( h ̂ KS ) | χ ), ξ f ( w ) ( r ) , and ρf(r).30 

  5. For each window and for each stochastic orbital, generate and store on the grid the stochastic projected orbitals ζ ( w ) ( r ) = r | ρ ̂ ( μ * ) P ̂ w | χ using the chemical potential determined in the previous step.

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

  7. Update the density and the KS Hamiltonian using the iterative subspace (DIIS) method37 and repeat the above steps until self-consistency is achieved, using the same random number seed.

We demonstrate the ew-efsDFT for a low G-center defect concentration in bulk silicon.38,39 We focus on the A-type38 G-center impurity embedded within a Si512 supercell (see the structure in Fig. 2). We performed Γ point DFT calculations using the Perdew–Burke–Ernzerhof (PBE)40 functional with Troullier–Martins norm-conserving pseudopotentials41 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 real-space grid spacing of 0.18 Å. In-gap states require a large β ≈ 900 in inverse Hartree to sufficiently converge the ground-state properties with respect to the electron temperature. 80 stochastic orbitals were used in both ew-efsDFT and o-efsDFT. The dressed fragments were Si64 with periodic boundary conditions, while the size of each core region was Si8. We used 41 energy windows in ew-efsDFT 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 orbitals26 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.

FIG. 2.

An A-type G-center (two carbon atoms shown in blue) embedded in a Si512 supercell (Si–Si bonds shown in yellow). The A-type G-center is constituted by a substitutional carbon atom, an interstitial carbon atom, and an interstitial silicon atom, which are highlighted as spheres.

FIG. 2.

An A-type G-center (two carbon atoms shown in blue) embedded in a Si512 supercell (Si–Si bonds shown in yellow). The A-type G-center is constituted by a substitutional carbon atom, an interstitial carbon atom, and an interstitial silicon atom, which are highlighted as spheres.

Close modal
FIG. 3.

SCF convergences from three single point calculations with different random number seeds are shown. In each single point calculation, the convergence is measured as the absolute difference between the energy per electron in the ith step (E) and the converged energy per electron (Ec). Units of energies are Hartree.

FIG. 3.

SCF convergences from three single point calculations with different random number seeds are shown. In each single point calculation, the convergence is measured as the absolute difference between the energy per electron in the ith step (E) and the converged energy per electron (Ec). Units of energies are Hartree.

Close modal

In Fig. 4, we assess the accuracy of ew-efsDFT for the electron density (upper panel) and the standard deviation in the electron density (lower panel) for selected positions in the vicinity of the G-center. The density and standard deviation were calculated from 5 independent ew-efsDFT or o-efsDFT 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 ew-efsDFT (shown in red) are significantly smaller than the corresponding o-efsDFT results (shown in blue). The noise of the electron density is significantly smaller, by approximately, a factor of 5 when compared to o-efsDFT.

FIG. 4.

Upper panel: the absolute value of the electron density difference (|Δρ|) between the stochastic (ew-efsDFT in red, o-efsDFT in blue) and a deterministic calculation. Lower panel: the standard deviation of electron density (ρSTD) evaluated with ew-efsDFT (red line) and o-efsDFT (blue line). |Δρ| and ρSTD are shown along the x axis with y = 5.3 Å and z = 5.0 Å. The peak marked by “Carbon” corresponds to a region closed to a carbon atom.

FIG. 4.

Upper panel: the absolute value of the electron density difference (|Δρ|) between the stochastic (ew-efsDFT in red, o-efsDFT in blue) and a deterministic calculation. Lower panel: the standard deviation of electron density (ρSTD) evaluated with ew-efsDFT (red line) and o-efsDFT (blue line). |Δρ| and ρSTD are shown along the x axis with y = 5.3 Å and z = 5.0 Å. The peak marked by “Carbon” corresponds to a region closed to a carbon atom.

Close modal

In Table I, we list the kinetic energy (Ek), the nonlocal pseudopotential energy (Enl), the Hartree energy (EH), the local pseudopotential energy (Eloc), the exchange–correlation energy (Exc), and the total energy (Etot), 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 ew-efsDFT compared to o-efsDFT 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 ew-efsDFT and o-efsDFT are slightly outside one standard deviation from the deterministic DFT result.30,31 We note in passing that the ew-sDFT approach (without fragmentation) does not reduce the noise in the total energy per electron, as discussed previously.32 

TABLE I.

The kinetic energy (Ek), the nonlocal pseudopotential energy (Enl), the Hartree energy (EH), the local pseudopotential energy (Eloc), the exchange–correlation energy (Exc), and the total energy (Etot), all per electron in eV obtained by deterministic DFT (dDFT), o-efsDFT, and ew-efsDFT. The standard deviation in the last digits of the energy is presented in parentheses.

Method Ek/Ne Enl/Ne EH/Ne Eloc Exc Etot/Ne
dDFT  10.1767  5.7871  2.0194  −8.7085  −8.1027  −26.8197 
o-efsDFT  10.1735(35)  5.7882(36)  2.0222(19)  −8.7074(40)  −8.1052(8)  −26.8205(11) 
ew-efsDFT  10.1778(11)  5.7874(4)  2.0193(6)  −8.7086(11)  −8.1028(3)  −26.8186(4) 
Method Ek/Ne Enl/Ne EH/Ne Eloc Exc Etot/Ne
dDFT  10.1767  5.7871  2.0194  −8.7085  −8.1027  −26.8197 
o-efsDFT  10.1735(35)  5.7882(36)  2.0222(19)  −8.7074(40)  −8.1052(8)  −26.8205(11) 
ew-efsDFT  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 x-direction for selected atoms obtained by the ew-efsDFT (upper panel) and the o-efsDFT (lower panel). Error bars indicate the standard deviation for each force. Clearly, the statistical fluctuations are significantly smaller for ew-efsDFT compared to o-efsDFT. In order to estimate the overall noise reduction efficiency, we averaged the standard deviations of Fx (nuclei force along the x axis) over all atoms. The averaged standard deviation of Fx is ≈0.53 eV/Å for o-efsDFT and ≈0.09 eV/Å for ew-efsDFT. Similar results were also obtained for other force components Fy and Fz, which implies that ew-efsDFT standard deviations of nuclei forces are about a factor of 6 smaller than those in o-efsDFT. In other words, to achieve a similar noise level in o-efsDFT would require ≈30 times more stochastic orbitals. No bias was observed for forces on the nuclei in ew-efsDFT.

FIG. 5.

Forces on the nuclei along the x-direction (Fx) for selected atoms calculated from ew-efsDFT (upper panel) and o-efsDFT (lower panel). Error bars in the forces on the nuclei were obtained from five runs. Blue symbols are Fx calculated from a deterministic DFT. The dashed line signifies the boundary between carbon and silicon atoms.

FIG. 5.

Forces on the nuclei along the x-direction (Fx) for selected atoms calculated from ew-efsDFT (upper panel) and o-efsDFT (lower panel). Error bars in the forces on the nuclei were obtained from five runs. Blue symbols are Fx calculated from a deterministic DFT. The dashed line signifies the boundary between carbon and silicon atoms.

Close modal

In Table II, we summarize the total wall time and the contribution from the Chebyshev moments, χ | T n ( h ̂ KS ) | χ , 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 ew-efsDFT and o-efsDFT methods. Each SCF iteration is ≈50% longer in ew-efsDFT, but the number of SCF iterations required to achieve a similar convergence is smaller in ew-efsDFT, resulting in similar wall times. Note that the statistical error in ew-efsDFT is much smaller than o-efsDFT. To achieve similar statistical errors in o-efsDFT, it would result in wall times that are roughly 30 times longer than ew-efsDFT.

TABLE II.

Averaged computational time (in hours) for o-efsDFT and ew-efsDFT. The wall time and number of SCF iterations are averaged over five independent runs. The time of calculating |ξ(w)⟩, |ζ(w)⟩, and |ξ⟩ and the time of projecting |ξ(w)⟩ with ρ ̂ frag are averaged over all SCF iterations and all five runs. All calculations were performed on a 40 node cluster computer, where each node contains two 16-core Intel Xeon Processors E5-2698 v3 at 2.3 GHz.

o-efsDFT ew-efsDFT
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 χ | T n ( h ̂ KS ) | χ (h)  0.32  0.73 
Time for generating |ξ⟩ or |ξ(w)⟩ (h)  0.64  0.72 
Time for projecting |ζ(w)⟩ with ρ ̂ frag (h)  N/A  0.04 
o-efsDFT ew-efsDFT
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 χ | T n ( h ̂ KS ) | χ (h)  0.32  0.73 
Time for generating |ξ⟩ or |ξ(w)⟩ (h)  0.64  0.72 
Time for projecting |ζ(w)⟩ with ρ ̂ frag (h)  N/A  0.04 

The main difference between the two methods is the computational time for generating the Chebyshev moments, χ | T n ( h ̂ KS ) | χ . In o-efsDFT, we used the relation T 2 n = 2 T n 2 1 to evaluate T n ( h ̂ KS ) | χ for n > Nc/2, thereby computing only 1/2 the number of moments. This relation cannot be used in ew-efsDFT 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 Nw projected stochastic orbitals in ew-efsDFT compared to only one projected orbital in o-efsDFT, resulting in ≈20% increase in generating |ξ⟩ vs |ξ(w)⟩. In addition, in ew-efsDFT, one has to compute | ξ f ( w ) with the fragment density matrix, i.e., i = 1 N occ f φ i f ( r ) φ i f | ξ ( w ) D f , in each SCF iteration, while in o-efsDFT, the projection of |χ⟩ is performed only once at the beginning of the calculation.

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 embedded-fragmented stochastic DFT27 with the energy window stochastic DFT.32 The new approach builds on both real-space and energy-space fragmentation, resulting in a significant reduction in the noise in single-particle observables without affecting the computational time. The performance of the ew-efsDFT was tested for a G-center embedded in bulk silicon with a small fundamental gap and in-gap impurity states, making this a rather challenging system for DFT. Compared to o-efsDFT and ew-DFT (not shown explicitly here), the statistical error in the forces is ∼6 times smaller in ew-efsDFT, 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 ew-efsDFT method to structural minimization is currently under way.33 

The asymptotic behavior of single-particle observables and analysis of the variance of the electron density for ew-efsDFT are provided in the supplementary material.

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. DE-AC02-05CH11231 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. DE-AC02-05CH11231. R.B. gratefully acknowledges support from the Germany-Israel Foundation (GIF) (Grant No. I-26-303.2-2018).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
P.
Hohenberg
and
W.
Kohn
,
Phys. Rev.
136
,
B864
(
1964
).
2.
W.
Kohn
and
L. J.
Sham
,
Phys. Rev.
140
,
A1133
(
1965
).
3.
A.
Aliano
and
G.
Cicero
, “
Ab initio DFT simulations of nanostructures
,” in
Encyclopedia of Nanotechnology
, edited by
B.
Bhushan
(
Springer Netherlands
,
Dordrecht
,
2012
), pp.
11
17
.
4.
T.
Frauenheim
,
G.
Seifert
,
M.
Elstner
,
T.
Niehaus
,
C.
Köhler
,
M.
Amkreutz
,
M.
Sternberg
,
Z.
Hajnal
,
A. D.
Carlo
, and
S.
Suhai
,
J. Phys.: Condens. Matter
14
,
3015
(
2002
).
5.
D. J.
Cole
and
N. D. M.
Hine
,
J. Phys.: Condens. Matter
28
,
393001
(
2016
).
7.
8.
T.
Zhu
,
W.
Pan
, and
W.
Yang
,
Phys. Rev. B
53
,
12713
(
1996
).
9.
W.
Yang
and
T. S.
Lee
,
J. Chem. Phys.
103
,
5674
(
1995
).
10.
A. W.
Götz
,
S. M.
Beyhan
, and
L.
Visscher
,
J. Chem. Theory Comput.
5
,
3161
(
2009
).
11.
T. A.
Wesolowski
,
S.
Shedge
, and
X.
Zhou
,
Chem. Rev.
115
,
5891
(
2015
).
12.
J. D.
Goodpaster
,
T. A.
Barnes
, and
T. F.
Miller
,
J. Chem. Phys.
134
,
164108
(
2011
).
13.
C.
Huang
and
E. A.
Carter
,
J. Chem. Phys.
135
,
194104
(
2011
).
14.
F.
Mauri
,
G.
Galli
, and
R.
Car
,
Phys. Rev. B
47
,
9973
(
1993
).
15.
P.
Ordejón
,
D. A.
Drabold
,
M. P.
Grumbach
, and
R. M.
Martin
,
Phys. Rev. B
48
,
14646
(
1993
).
16.
S.
Goedecker
,
J. Comput. Phys.
118
,
261
(
1995
).
17.
E.
Hernández
and
M. J.
Gillan
,
Phys. Rev. B
51
,
10157
(
1995
).
19.
R.
Baer
and
M.
Head-Gordon
,
Phys. Rev. Lett.
79
,
3962
(
1997
).
20.
A. H. R.
Palser
and
D. E.
Manolopoulos
,
Phys. Rev. B
58
,
12704
(
1998
).
21.
J.
VandeVondele
,
U.
Borštnik
, and
J.
Hutter
,
J. Chem. Theory Comput.
8
,
3565
(
2012
).
22.
A.
Nakata
,
Y.
Futamura
,
T.
Sakurai
,
D. R.
Bowler
, and
T.
Miyazaki
,
J. Chem. Theory Comput.
13
,
4146
(
2017
).
23.
L. E.
Ratcliff
,
G. J.
Conduit
,
N. D. M.
Hine
, and
P. D.
Haynes
,
Phys. Rev. B
98
,
125123
(
2018
).
24.
T. D.
Kühne
,
M.
Iannuzzi
,
M.
Del Ben
,
V. V.
Rybkin
,
P.
Seewald
,
F.
Stein
,
T.
Laino
,
R. Z.
Khaliullin
,
O.
Schütt
,
F.
Schiffmann
,
D.
Golze
,
J.
Wilhelm
,
S.
Chulkov
,
M. H.
Bani-Hashemian
,
V.
Weber
,
U.
Borštnik
,
M.
Taillefumier
,
A. S.
Jakobovits
,
A.
Lazzaro
,
H.
Pabst
,
T.
Müller
,
R.
Schade
,
M.
Guidon
,
S.
Andermatt
,
N.
Holmberg
,
G. K.
Schenter
,
A.
Hehn
,
A.
Bussy
,
F.
Belleflamme
,
G.
Tabacchi
,
A.
Glöß
,
M.
Lass
,
I.
Bethune
,
C. J.
Mundy
,
C.
Plessl
,
M.
Watkins
,
J.
VandeVondele
,
M.
Krack
, and
J.
Hutter
,
J. Chem. Phys.
152
,
194103
(
2020
).
25.
L.
Lin
,
J.
Lu
,
L.
Ying
, and
W.
E
,
J. Comput. Phys.
231
,
2140
(
2012
).
26.
R.
Baer
,
D.
Neuhauser
, and
E.
Rabani
,
Phys. Rev. Lett.
111
,
106402
(
2013
).
27.
M.
Chen
,
R.
Baer
,
D.
Neuhauser
, and
E.
Rabani
,
J. Chem. Phys.
150
,
034106
(
2019
).
28.
D.
Neuhauser
,
R.
Baer
, and
E.
Rabani
,
J. Chem. Phys.
141
,
041102
(
2014
).
29.
E.
Arnon
,
E.
Rabani
,
D.
Neuhauser
, and
R.
Baer
,
J. Chem. Phys.
146
,
224111
(
2017
).
30.
M. D.
Fabian
,
B.
Shpiro
,
E.
Rabani
,
D.
Neuhauser
, and
R.
Baer
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
9
,
e1412
(
2019
).
31.
Y.
Cytter
,
E.
Rabani
,
D.
Neuhauser
, and
R.
Baer
,
Phys. Rev. B
97
,
115207
(
2018
).
32.
M.
Chen
,
R.
Baer
,
D.
Neuhauser
, and
E.
Rabani
,
J. Chem. Phys.
151
,
114116
(
2019
).
33.
M.
Chen
,
R.
Baer
,
D.
Neuhauser
, and
E.
Rabani
, “
Geometry Minimization with Stochastic Density Functional Theory
” (unpublished).
34.
R.
Kosloff
,
J. Phys. Chem.
92
,
2087
(
1988
).
35.
R.
Kosloff
,
Annu. Rev. Phys. Chem.
45
,
145
(
1994
).
36.
F.
Shimojo
,
R. K.
Kalia
,
A.
Nakano
, and
P.
Vashishta
,
Phys. Rev. B
77
,
085103
(
2008
).
37.
38.
L. W.
Song
,
X. D.
Zhan
,
B. W.
Benson
, and
G. D.
Watkins
,
Phys. Rev. B
42
,
5765
(
1990
).
39.
E.
Rotem
,
J. M.
Shainline
, and
J. M.
Xu
,
Appl. Phys. Lett.
91
,
051127
(
2007
).
40.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
41.
N.
Troullier
and
J. L.
Martins
,
Phys. Rev. B
43
,
1993
(
1991
).
42.
L.
Kleinman
and
D. M.
Bylander
,
Phys. Rev. Lett.
48
,
1425
(
1982
).

Supplementary Material