Over this past decade, we combined the idea of stochastic resolution of identity with a variety of electronic structure methods. In our stochastic Kohn-Sham density functional theory (DFT) method, the density is an average over multiple stochastic samples, with stochastic errors that decrease as the inverse square root of the number of sampling orbitals. Here, we develop a stochastic embedding density functional theory method (se-DFT) that selectively reduces the stochastic error (specifically on the forces) for a selected subsystem(s). The motivation, similar to that of other quantum embedding methods, is that for many systems of practical interest, the properties are often determined by only a small subsystem. In stochastic embedding DFT, two sets of orbitals are used: a deterministic one associated with the embedded subspace and the rest, which is described by a stochastic set. The method agrees exactly with deterministic calculations in the limit of a large number of stochastic samples. We apply se-DFT to study a p-nitroaniline molecule in water, where the statistical errors in the forces on the system (the p-nitroaniline molecule) are reduced by an order of magnitude compared with nonembedding stochastic DFT.

DFT (Density Functional Theory) traditionally follows the Kohn-Sham scheme where a set of one-particle equations is solved self-consistently. For large systems, the solution of these equations scales as O(Ne2)O(Ne3) with the number of electrons Ne, so there is a lot of interest in variants that scale linearly with system size. Such methods include orbital-free DFT with density-dependent kinetic energy functionals,1,2 linear-scaling approaches where the system is split into parts that are woven together via constraints,3 and embedding techniques where an inner part is treated by DFT and an outer part by orbital-free DFT.4–6 

Previously, we developed stochastic DFT (sDFT), a method that can be viewed as a bridge between Kohn-Sham DFT and orbital-free DFT.7,8 Instead of computing Kohn-Sham orbitals for all occupied states, we apply a Chebyshev filter to a few stochastic orbitals7 and extract the density from these filtered orbitals, circumventing the time-consuming diagonalization step. In the limit of infinitely many stochastic samples, the stochastic errors approach zero and the results agree exactly with deterministic calculations.

In practice, the sample size will be finite, and there will be stochastic noises associated with the results. There is no problem with noisy results per se. People use noise to simulate temperature effects. For example, the noisy forces of sDFT have been used to determine the structure of nanocrystals at finite temperatures by sampling the Boltzmann distribution through the technique known as Langevin dynamics.9 When the noisy forces have a large standard deviation, the friction coefficient in the Langevin equation must be large and this deteriorates the sampling efficiency. Thus, there is a great advantage in reducing the fluctuations in the noise. Hence, our aim in the manuscript is to discover a way to reduce the noise in the force by a significant factor than in our previous applications.

In a follow-up work,9–11 we have shown how to reduce the standard deviation in sDFT (and therefore accelerate the convergence), using a method we label stochastic-fragment DFT (sf-DFT). Here, instead of sampling stochastically the full density, we sample stochastically only the difference between the full density and a zero-order density which is easy to calculate. The difference is generally small, thereby reducing the fluctuations.

Here, we develop an alternate method whereby a given subregion is embedded. Essentially, this subregion is treated deterministically, while the rest of the system is treated stochastically (this is a simplifying view, and the more precise methodology is described later). The motivation for this method is that, for many realistic systems, only a subsystem is of particular importance. The idea of embedding was widely adopted to treat such systems. In most embedding methods, the subsystem of interest is calculated at higher level theories, while the rest is treated with less accurate but more efficient methods to reduce computational cost.4,5,12–23 Our stochastic density functional theory embedding method adopts an analogous strategy, except that here the larger stochastic region embeds the smaller deterministic part.

An attractive feature of the stochastic embedding method is that the errors due to the embedding are numerically controlled since in the limit of infinitely many stochastic samplings, the method is exact. As such, there is no residual arbitrariness due to the choice of an embedding potential.

This paper is organized as follows. Sec. I presents the theory, and the practical algorithm is reviewed in Sec. II. In Sec. III, a practical system is studied, embedding of a dye (p-nitroaniline) in 216 water molecules. Discussion and possible extensions follow in Sec. IV.

We first review stochastic DFT as developed in our previous works.7,10–12

In DFT, the key component is the electron density ρr, which we express as the trace of a Heaviside step function,

ρr2=r|Θ(μH)|r,
(1)

where we assume spin-unpolarized DFT. Here, μ is the electron chemical potential, determined by ensuring the correct total number of electrons,

Ne=ρ(r)dr,
(2)

and the one-body Hamiltonian is H=122+v(r), where we introduced the effective one-electron potential due to the nuclear (vN) electron-electron Coulomb interaction (vH) and exchange-correlation (vXC) parts. We assume for simplicity that the exchange-correlation (and therefore the total effective) potential depends on the local density, v=vρ.

In the usual deterministic formulations of Kohn-Sham DFT, the electron density is expressed as the sum over one-electron states, and the total number of electrons is determined by the occupation number of each state. Thus, the Heaviside filter becomes a projection to the occupied subspace,

Θ(μH)=iNocc|ψiψi|,
(3)

where we introduced the number of occupied orbitals (Nocc = Ne/2). The one-electron orbitals ψi are obtained by diagonalization of the effective one-electron Hamiltonian matrix, resulting in a nominal Ne3 scaling of Kohn-Sham DFT. Expectation values of one-electron operators are obtained from the occupied states,

A=iNoccψi|A|ψi.
(4)

In stochastic Kohn-Sham DFT, on the other hand, we use a set of stochastic orbitals, ξ(r), with the property that

|ξξ|ξ=I,
(5)

where the curly brackets stand for averaging over all stochastic orbitals. Inserting the identity operator in the expression for ρ(r), the electron density is thus expressible as

ρ(r)=r|Θ12|ξξ|Θ12|rξ=|ξμ(r)|2ξ,
(6)

where we abbreviate Θ12Θ12(μH) and ξμΘ12ξ.

The filtered stochastic orbitals are linear combinations of all occupied states with random coefficients,

ξμ(r)=iNoccciϕi(r),
(7)

where ci = ⟨ϕi|ξ⟩, so

ci*cj=δij.
(8)

Similarly, the expectation values of any one-particle operator D is

D=ξμ|D|ξμξ1Nsξξμ|D|ξμ.
(9)

Here, Ns is the number of stochastic orbitals used in practice. As in any stochastic method, the expectation value, obtained as the average of a finite number of samples, will have an associated stochastic error which is proportional to 1/Ns. The actual number of stochastic orbitals is chosen based on the required level of precision.

In practice, the method relies on the fact that the calculation of each stochastic vector scales only linearly with system size. Specifically, we use a smooth Heaviside function, Θ(μH)=12erfcβμH, where β needs to be much larger than the inverse band gap. The smooth theta function is then expressed as a finite sum of Chebyshev polynomials, Θ12=nan(μ)Tn(Hscaled), where Hscaled is a scaled Hamiltonian with eigenvalues in the range 1,1 and an(μ) are the Chebyshev coefficients of 12erfcβμH12. Therefore,

|ξμ=nan(μ)|ξn,
(10)

where the Chebyshev vectors are obtained recursively, |ξn⟩ = 2Hscaled|ξn−1⟩ − |ξn−2⟩, and |ξn=0⟩ ≡ |ξ⟩.

The Chebyshev expansion makes it possible to analytically determine the chemical potential. Specifically, expand Θ = ∑nbn(μ)Tn(Hscaled), where bn(μ) are the Chebyshev coefficients of 12erfcβμH. Then, using Eq. (2) gives

Ne2=1Nsξξμ|Θ(μH)|ξμ=nbn(μ)Rn,
(11)

where Rn=Ns1ξξ|TnHscaled|ξ. Therefore, μ is varied until Eq. (11) is fulfilled.

Next, we turn to embedding, first traditional and then stochastic.

As in other embedding methods, the motivation for stochastic DFT embedding is that in many practical applications, the properties of a system depend on much smaller subsystem(s), such as defects in semiconductors or active sites of proteins. Often, we do not even care for the rest of the system except the embedded part. Even when a quantum-mechanical treatment of the rest of the system (the environment) is necessary, the level of precision required is usually not as high as that of the subsystem(s) of interest.

A key component in all types of embedding methods is the specific quantity through which the properties of the environment are conveyed to the subsystem(s). In DFT embedding, the quantity is the electron density of the entire system.18 Specifically, subset A denotes the subsystem of interest, and subset B is associated with the environment (the precise meaning of these two subsets would be flexible). The total electron density is therefore partitioned ρ = ρA + ρB.

In traditional deterministic DFT embedding, one then derives an approximate functional for the two regions that captures the energy of the environment as well as its interaction with the subsystem(s) of interest. Solving this equation for orbitals in subset A is equivalent to solving an ordinary Kohn-Sham equation with an extra external potential due to the embedding. Following Kohn Sham formulation, the external potential can be written as

vA=δδρA(Ts[ρ]Ts[ρA])+VH[ρB]+δδρA(Exc[ρ]Exc[ρA]).

Of all the three terms to the right, the first term is referred to as the nonadditive kinetic potential, which usually has the dominating contribution. Heuristically, the term can be thought of arising from the requirement that the orbitals of the subsystems(s) should be orthogonal to the orbitals of the environment.

In practice, this term can either be derived from its functional form,4,18 or the orthogonality can be imposed.12–15,24,25

In our stochastic embedding DFT (se-DFT), there is the same loose overall goal as in deterministic embedding, i.e., the different treatment of a smaller system and the environment. However, the methods are quite different. In our approach, the embedded space is treated deterministically and the other (environment) is treated stochastically but at the same level of computational methods. Therefore, the only sense in which embedding is approximate here is numerical, i.e., if we use enough stochastic orbitals, the results will agree exactly with fully deterministic calculations. The two spaces (system and environment) see the same overall Hamiltonian, and there is no uncontrolled ansatz.

Specifically, using the same language of partitioning space to parts, we separate the total electron density into two parts, abbreviating

ρ(r)2=r|Θ12Θ12|r=r|Θ12PΘ12|r+r|Θ12QΘ12rρA(r)2+ρB(r)2.
(12)

The first term in the splitting projects onto the “A” subspace, defined by its basis

P^=iA|χiχi|,       χi|χj=δij.
(13)

Therefore,

12ρA(r)=iAr|Θ12|χiχi|Θ12|r=iA|χi,μ(r)|2,
(14)

where χi,μ=Θ12(μH)χi. Note that the χi(r) basis does not have to be related directly to the molecular orbitals but would typically be from a set of atomic orbitals in a given region although there is a lot of freedom in the definition. For example, in the example studied later, we choose a set of local Gaussian atomic orbitals on each atom in the embedded subsystem [labeled ϕi(r), iA] and then orthogonalize them, to produce χi(r)=j(S12)ijϕj(r), where Sij = ⟨ϕi|ϕj⟩ is the overlap matrix of the embedded-part atomic orbitals.

The second term in the splitting is associated with QIP, the orthogonal projection to the other (“B”) subspace. Since Q2 = Q and inserting the identity operator I=|ξξ|ξ, we get

12ρB(r)=r|Θ12QQΘ12|r=r|Θ12Q|ξξ|QΘ12|rξ=|ξ¯μ(r)|2ξ,
(15)

where

|ξ¯μΘ12μHQ|ξ
(16)

is obtained by two consecutive projections: first a random orbital is projected to the space orthogonal to the embedded P part, and the result is then projected to the occupied space of the full system.

Thus, we reach the main embedding expression, the separation of the density into two parts,

12ρ(r)=iA|χi,μ(r)|2+1Nsξ|ξ¯μ(r)|2,
(17)

one associated with the deterministic subspace and one with an orthogonal stochastic part. The two parts are connected through the application of the density matrix operator, ΘμH, since the potential in the Hamiltonian depends on the density, v=v[ρ], and the density is a mixture of stochastic and deterministic parts.

An important feature of the algorithm is that the deterministic and stochastic orbitals, χ and ξ¯μ, that make up the density [Eq. (17)] are not orthogonal—neither among themselves nor to each other. The orthogonality of the P and Q spaces reflects in the orthogonality of the original χi and functions, but that orthogonality is lost when we act on χi and on with Θ12(μH) in Eqs. (14) and (16). Further note that, as mentioned, the same overall Hamiltonian [and therefore the same Θ12(μH)] is used in preparing both the stochastic and deterministic orbitals, i.e., they are treated on equal footing.

The procedure described above serves to introduce a general and novel paradigm of embedding deterministic and stochastic descriptions together. The paradigm is not limited to DFT and would be generally applicable to other techniques, for example, time-dependent density functional theory (TDDFT), so that the inner region would be described by deterministically propagated orbitals, or other methods.

We also note that the approximation in this approach, i.e., the stochastic part, is numerically controllable. With more stochastic orbitals, the error decreases. It does present a fully quantum alternative to other embedding methods, where one needs to increase the size of the embedded part to ensure convergence. Furthermore, the method scales well for a large method. The main cost is proportional to NstochNembed2, so it is very manageable even for very large regions with thousands of basis set functions in the embedded region. Like all stochastic methods, the aim is giant systems, where the gentle scaling of the method enables a calculation.

The overall stochastic embedding DFT method is then quite similar to the stochastic DFT algorithm:

  1. Generate Ns stochastic orbitals: ξ(r)=±1/d3r, where d3r is the volume element associated with the grid. Also create a reasonable initial density ρ(r), which integrates to the correct number of valence electrons.

  2. Determine the one-particle effective potential and Hamiltonian H = T + v[ρ].

  3. For each stochastic orbital, project out the components along the atomic basis functions, i.e., prepare ξ¯=Qξ,

ξ¯(r)=ξ(r)iciχi(k),

where

ci=χi(r)ξ(r)dr.
  1. Determine the correct chemical potential μ as the one that integrates correctly the total density, i.e., from Eq. (11) where now

Rn=iAχi|TnHscaled|χi+1Nsξξ¯|TnHscaled|ξ¯,
(18)

i.e., the residues and therefore the constraint on the integrated density include both the deterministic and stochastic parts.

  1. Chebyshev filter the orthogonalized atomic basis functions as well as the projected stochastic functions,

|χi,μ=Θ12(μH)|χi,|ξ¯μ=Θ12(μH)|ξ¯,
  1. Calculate the charge densities for this μ from Eq. (17).

  2. Reiterate steps 2–6 until the density does not vary, i.e., SCF convergence is reached.

With the filtered atomic basis functions and filtered stochastic orbitals, the expectation value of any one-particle operator D is

D=iχi,μ|D|χi,μ+1Nsξ¯μ|D|ξ¯μ.
(19)

The algorithm is therefore very similar to the original stochastic DFT approach. The only differences are that (i) in addition to the stochastic functions, one also needs to project the density matrix [or more precisely Θ12μH] on the deterministic basis making the embedded part, and (ii) for the stochastic part, we now project out the deterministic part (i.e., apply Q) before filtering with the Chebyshev expansion of Θ12μH.

We applied the stochastic density functional embedding method to study a realistic case of embedding, i.e., a dye in water. The dye was a p-nitroaniline molecule, and it was embedded in 216 water molecules.

The configuration of the system was obtained from snapshots of molecular dynamic (MD) simulations with Gromacs 5.26 The dynamics simulations used a generalized amber force field with charges from AM1-BCC for p-nitroaniline27 and TIP4P with allowed flexibility of bond/bend for water.28 

The simulation steps for the preparation of the configuration were standard, involving first high temperature equilibrations of the p-nitroaniline, followed by NVT simulations at room temperature, and then NPT equilibration at room temperature and pressure. We then ran the MD simulation and sampled a specific configuration after several nanoseconds. This configuration was used for subsequent DFT calculations and is shown in Fig. 1.

FIG. 1.

Configuration of the p-nitroaniline/water system used in DFT calculation. The p-nitroaniline molecule is represented by balls and sticks, and water molecules are represented by wire-frames.

FIG. 1.

Configuration of the p-nitroaniline/water system used in DFT calculation. The p-nitroaniline molecule is represented by balls and sticks, and water molecules are represented by wire-frames.

Close modal

We note that in a comprehensive study of the solvent effect, the molecular properties should be averaged over multiple configurations. However, as the purpose of the current study is to demonstrate the capacity of the stochastic embedding DFT method, we used single configuration instead. The restriction also allows us to focus on the difference in results obtained from the embedding method and purely stochastic method.

For the DFT calculation, we imposed periodic boundary conditions. A plane wave expansion was used, with atomic norm-conserving Troullier-Martins pseudopotentials replacing the core-valence interaction. A local-density approximation (LDA) functional was used. An 883 grid with a spacing of 0.402 atomic units was used, while the plane wave kinetic-energy cutoff was 15 hartree. The inverse-temperature-like parameter β in the Heaviside function erfcβμH was set at β = 0.03 hartree−1, requiring 1173 Chebyshev propagations in acting with Θ12(μH).

For the basis functions {|χi⟩}, we used a Gaussian double-zeta basis optimized for pseudopotentials, as given in the Quickstep29 data set.30 

Three sets of calculations were carried out. The first used Ns = 96 stochastic orbitals, without embedding. The second used the same Ns = 96 stochastic functions but supplemented them with the double zeta atomic basis set for all 16 atoms belonging to the p-nitroaniline molecule. This deterministic basis set for the dye contained 160 functions (an average of 10 functions per atom). The atomic functions were orthogonalized, giving rise to 160 orthogonal χi(r) functions. Therefore, a total of 256 functions were employed (160 deterministic and 96 stochastic).

Since the number of deterministic functions in the second, embedded, set of calculations was quite large, we also compared the second set with a third set where all functions were stochastic, and where 256 functions were used, i.e., the same overall number as in the second set. Thus, the numerical effort, mostly associated with the Chebyshev application of Θ12, is similar in the second and third sets.

Each set of calculations was repeated ten times, with different random seeds for generating the stochastic orbitals, to obtain the standard deviation of the stochastic approaches.

To benchmark our results, we also performed a conventional deterministic Kohn-Sham DFT calculation which matches the results of Quantum-Espresso.31 

The most time-consuming step in the stochastic DFT formulation is the application of the Chebyshev filter. Therefore, as mentioned, the time required to perform calculation with embedding and Ns = 96 (second set) is comparable to that with no embedding and Ns = 256 (third set) and is about 2.5 times the time required to perform calculation with no embedding and Ns = 96 (first set). Indeed, in practice, about 14 core hours per SCF iteration (on a cluster with 2.5 GHz nodes) were needed for the second and third sets, with about 6 core hours for the first set. The deterministic set required about 11 core hours per iterations. In all cases, 30 direct inversion of the iterative subspace iterations were used for full SCF convergence.

Next, we compared the total energy per electron obtained from the three sets of calculations with that obtained from the deterministic calculation, as well as the individual contribution from Hartree and exchange-correlation energies. The comparison is shown in Table I. The results show good overall agreement between the stochastic and deterministic calculations. In the table, we also included standard deviations of the various results. The first observation is that results obtained from stochastic embedding calculations are consistent with purely stochastic calculations of the same NS. This observation is expected because the various terms contributing to total energies are dominated by the 216 water molecules. Another observation is that the standard deviation can be made smaller by increasing NS. In fact, the standard deviation is expected to be proportional to 1/NS.

TABLE I.

DFT energies per electron, in Hartree. For the stochastic calculations, the energies were obtained as average of ten calculations, and the standard deviation is included. As seen in the table, results obtained from stochastic embedding calculations agree well with stochastic calculations of the same NS.

Ns = 96 withoutNs = 96 withNs = 256 without
embedding embeddingembeddingDeterministic
Total energy per electron −2.115 ± 0.0015 −2.115 ± 0.0016 −2.117 ± 0.0008 −2.117 
Hartree energy per electron 1.098 ± 0.0023 1.098 ± 0.0025 1.100 ± 0.0016 1.103 
Exchange-correlation energy per electron −0.521 ± 0.0004 −0.521 ± 0.0004 −0.521 ± 0.0002 −0.522 
Ns = 96 withoutNs = 96 withNs = 256 without
embedding embeddingembeddingDeterministic
Total energy per electron −2.115 ± 0.0015 −2.115 ± 0.0016 −2.117 ± 0.0008 −2.117 
Hartree energy per electron 1.098 ± 0.0023 1.098 ± 0.0025 1.100 ± 0.0016 1.103 
Exchange-correlation energy per electron −0.521 ± 0.0004 −0.521 ± 0.0004 −0.521 ± 0.0002 −0.522 

The effect of embedding comes to play in quantities relating to the embedded region, and such a main quantity is the force on each atom. Figure 2 shows the overall forces and their standard deviations for 60 out of the 644 atoms in the sample; the first 16 are the embedded dye, and the other 44 are from the water molecules. The figure clearly shows that, with embedding, the forces on the embedded atoms have much higher accuracy (much smaller standard deviation) than the forces on the water molecule.

FIG. 2.

Forces and their standard deviations for 60 (out of the 664 overall) atoms, in eV/Å. The first 16 atoms in the plot are the p-nitroaniline molecule, and the other 44 are from water molecules. In the second set, the atoms of the p-nitroaniline are embedded, and in the first and third set, they are not. The blue circles denote the forces calculated by deterministic DFT, while the center of each bar refers to the average stochastic force. Red error bars are associated with the atoms of the p-nitroaniline molecule, and green bars are used for the standard deviation of the forces on the water atoms. Plot a) shows the results of using 96 stochastic orbitals without embedding; plot b) shows the results of using 96 stochastic orbitals, plus 160 basis functions used for the embedded p-nitroaniline molecule; while plot c) shows the results of using 256 stochastic orbitals, without embedding.

FIG. 2.

Forces and their standard deviations for 60 (out of the 664 overall) atoms, in eV/Å. The first 16 atoms in the plot are the p-nitroaniline molecule, and the other 44 are from water molecules. In the second set, the atoms of the p-nitroaniline are embedded, and in the first and third set, they are not. The blue circles denote the forces calculated by deterministic DFT, while the center of each bar refers to the average stochastic force. Red error bars are associated with the atoms of the p-nitroaniline molecule, and green bars are used for the standard deviation of the forces on the water atoms. Plot a) shows the results of using 96 stochastic orbitals without embedding; plot b) shows the results of using 96 stochastic orbitals, plus 160 basis functions used for the embedded p-nitroaniline molecule; while plot c) shows the results of using 256 stochastic orbitals, without embedding.

Close modal

As a side remark note that the forces on the nonembedded atoms have the same overall statistical fluctuations as without embedding. We know from previous studies that for large systems such as the present one, the stochastic errors are independent of size and depend only on the number of stochastic samples. Thus, the approach presented here would not deteriorate with the overall size of the full system.

Coming back to the embedded system (the dye), the higher accuracy on the dye forces is shown more quantitatively in Table II, where the standard deviations of the forces on the dye are an order of magnitude smaller than for water molecules.

TABLE II.

Average standard deviations of the atomic force magnitudes, in eV/Å.

Ns = 96 withoutNs = 96 withNs = 256 without
embedding embedding embedding
Average over 1.92 0.20 1.54 
embedded atoms    
Average over 2.23 2.22 1.47 
other atoms    
Ns = 96 withoutNs = 96 withNs = 256 without
embedding embedding embedding
Average over 1.92 0.20 1.54 
embedded atoms    
Average over 2.23 2.22 1.47 
other atoms    

The results show that embedding significantly reduces the stochastic error of the accelerations for the selected (i.e., embedded) atoms. As expected, when the same number of stochastic orbitals is used, the standard deviation averaged over the nonembedded atoms remains the same. Meanwhile, due to the good description of the embedded atoms by the deterministic atomic basis, the standard deviations of the forces for those atoms decrease by one order of magnitude relative to the no-embedding case. To achieve without embedding the same level of accuracy in the forces, we will need 10 000 stochastic orbitals, which would be very time-consuming.

In summary, we presented a stochastic embedding DFT method (se-DFT) that significantly reduces the statistical errors in the forces for the selected subgroup of atoms (i.e., the embedded atoms). Combined with the favorable linear scaling of stochastic DFT, the method can be applied to large systems of practical interests. Of course, as it stands, the method is not efficient for overall MD of the full system, due to the large stochastic errors on the environment (i.e., nonembedded) atoms; rather, it is suitable for applications where information on a selected region is desired. Results showed that the stochastic embedding method can indeed selectively reduce the stochastic errors associated with quantities on a specific subsystem. In the current work, we focused on the forces on embedded atoms, and we expect the same to hold for other quantities as well.

The embedding approach presented here is very general and can be extended in several directions. First, here we used an embedded space made from low-level atomic basis functions; we can replace it by a more general higher level basis. Furthermore, we could economize and choose in the P basis only occupied eigenfunctions of the embedded part in a dielectric medium. There will be occasions where the best basis would be energy selective, i.e., a few energy-selective molecular eigenfunctions or a few energy selective eigenfunctions from a large cluster would be best used.

A second direction is a combination of embedding with our previous overlapping fragment technique (sf-DFT). This method reduces the statistical error of stochastic DFT calculations9,10 for all atoms, typically by up to an order of magnitude; specifically, instead of stochastically sampling the full density, we sample the difference between the full density and a zeroth-order density, ρ(r), which is a solution of a simple zeroth-order Hamiltonian H0 (e.g., that of overlapping fragments). Specifically,

ρ(r)=ρ0(r)+|r|Θ12|ξ|2|r|Θ012|ξ|2ξ,
(20)

where Θ0 ≡ Θ(μ0H0) and μ0 is arbitrary. It is clear that this overlapping-fragment definition can be further extended by inserting projection operators as done earlier in the paper for the original stochastic DFT method. In a future paper, we will examine whether a combination of overlapping fragments with embedding (i.e., combining se-DFT with sf-DFT) reduces the errors in the forces of the embedded part even further than either method alone.

A third direction is for methods other than DFT. For example, embedding is applicable for the sublinear scaling stochastic TDDFT method developed by us.32 It is straightforward to see that the main embedding equation, Eq. (17), follows straightforwardly to TDDFT except that now all quantities (the density and the deterministic and stochastic orbitals) are time dependent. Such a time-dependent method has the desired property that the same Hamiltonian guides both the deterministic and stochastic functions. An embedding-TDDFT method would be applicable to study the change in optical properties of chromophores due to the presence of solvent molecules, where we can use only a few stochastic orbitals to sample the solvent molecules, while the chromophore will be treated with embedding. This direction would be explored in a future paper.

D.N. acknowledges support from the NSF, Grant No. CHE-1763176. E.R. acknowledges support from the Department of Energy, Photonics at Thermodynamic Limits Energy Frontier Research Center, under Grant No. DE-SC0019140. R.B. acknowledges support from the US-Israel Binational Science Foundation under the BSF-NSF program, Grant No. 2015687.

The work also used 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. The calculations were performed as part of the XSEDE33 computational Project No. TG-CHE170058.

Here, we give a simple demonstration of why embedding should improve the statistics. Say that the overall problem has two spaces that are essentially separate, so each eigenstate ϕi belongs to either the A or B subspaces. We choose the P subspace spanned by {χi : iA} such that it is close to the subspace spanned by states in A,

P|ϕicio|ϕi,
(A1)

where

cio1,iA,0,iB.
(A2)

The projected filtered stochastic orbital ξ¯μ(r) is then given by a linear combination,

ξ¯μ(r)=iNOcc(cicio)ϕi(r).
(A3)

For states belonging to A, instead of sampling ci, we are sampling (cicio) with average (1cio)0 and a much smaller standard deviation.

1.
W. C.
Witt
,
B. G.
del Rio
,
J. M.
Dieterich
, and
E. A.
Carter
, “
Orbital-free density functional theory for materials research
,”
J. Mater. Res.
33
(
07
),
777
795
(
2018
).
2.
J.
Lehtomäki
,
I.
Makkonen
,
M. A.
Caro
,
A.
Harju
, and
O.
Lopez-Acevedo
, “
Orbital-free density functional theory implementation with the projector augmented-wave method
,”
J. Chem. Phys.
141
(
23
),
234102
(
2014
).
3.
T.-S.
Lee
,
J. P.
Lewis
, and
W.
Yang
, “
Linear-scaling quantum mechanical calculations of biological molecules: The divide-and-conquer approach
,”
Comput. Mater. Sci.
12
(
3
),
259
277
(
1998
).
4.
T. A.
Wesolowski
and
A.
Warshel
, “
Frozen density functional approach for ab initio calculations of solvated molecules
,”
J. Phys. Chem.
97
(
30
),
8050
8053
(
1993
).
5.
T. A.
Wesolowski
,
S.
Shedge
, and
X.
Zhou
, “
Frozen-density embedding strategy for multilevel simulations of electronic structure
,”
Chem. Rev.
115
(
12
),
5891
5928
(
2015
).
6.
M.
Iannuzzi
,
B.
Kirchner
, and
J.
Hutter
, “
Density functional embedding for molecular systems
,”
Chem. Phys. Lett.
421
(
1-3
),
16
20
(
2006
).
7.
R.
Baer
,
D.
Neuhauser
, and
E.
Rabani
, “
Self-averaging stochastic Kohn-Sham density-functional theory
,”
Phys. Rev. Lett.
111
(
10
),
106402
(
2013
).
8.
Y.
Cytter
,
E.
Rabani
,
D.
Neuhauser
, and
R.
Baer
, “
Stochastic density functional theory at finite temperatures
,”
Phys. Rev. B
97
,
115207
(
2018
).
9.
E.
Arnon
,
E.
Rabani
,
D.
Neuhauser
, and
R.
Baer
, “
Equilibrium configurations of large nanostructures using the embedded saturated-fragments stochastic density functional theory
,”
J. Chem. Phys.
146
(
22
),
224111
(
2017
).
10.
D.
Neuhauser
,
R.
Baer
, and
E.
Rabani
, “
Communication: Embedded Fragment Stochastic Density Functional Theory
,”
J. Chem. Phys.
141
,
041102
(
2014
).
11.
M.
Chen
,
R.
Baer
,
D.
Neuhauser
, and
E.
Rabani
, “
Overlapped embedded fragment stochastic density functional theory for covalently-bonded materials
,”
J. Chem. Phys.
150
(
3
),
034106
(
2019
).
12.
F. R.
Manby
,
M.
Stella
,
J. D.
Goodpaster
, and
T. F.
Miller
, “
A simple, exact density-functional-theory embedding scheme
,”
J. Chem. Theory Comput.
8
(
8
),
2564
2568
(
2012
).
13.
P. K.
Tamukong
,
Y. G.
Khait
, and
M. R.
Hoffmann
, “
Density differences in embedding theory with external orbital orthogonality
,”
J. Phys. Chem. A
118
(
39
),
9182
9200
(
2014
).
14.
B.
Hégely
,
P. R.
Nagy
,
G. G.
Ferenczy
, and
M.
Kállay
, “
Exact density functional and wave function embedding schemes based on orbital localization
,”
J. Chem. Phys.
145
(
6
),
064107
(
2016
).
15.
C.
Tanner
,
K. R.
Brorsen
, and
S.
Hammes-Schiffer
, “
Communication: Density functional theory embedding with the orthogonality constrained basis set expansion procedure
,”
J. Chem. Phys.
146
(
21
),
211101
(
2017
).
16.
A. S.
Pereira Gomes
,
C. R.
Jacob
, and
L.
Visscher
, “
Calculation of local excitations in large systems by embedding wave-function theory in density-functional theory
,”
Phys. Chem. Chem. Phys.
10
(
35
),
5353
5362
(
2008
).
17.
D.
Loco
,
L.
Lagardère
,
S.
Caprasecca
,
F.
Lipparini
,
B.
Mennucci
, and
J.-P.
Piquemal
, “
Hybrid QM/MM molecular dynamics with AMOEBA polarizable embedding
,”
J. Chem. Theory Comput.
13
(
9
),
4025
4033
(
2017
).
18.
Q.
Sun
and
G. K.-L.
Chan
, “
Quantum embedding theories
,”
Acc. Chem. Res.
49
(
12
),
2705
2712
(
2016
).
19.
T. A.
Wesolowski
and
A.
Savin
, “
Non-additive kinetic energy and potential in analytically solvable systems and their approximated counterparts
,” in
Recent Progress in Orbital-free Density Functional Theory
(
World Scientific
,
2013
), pp.
275
295
.
20.
D. G.
Jason
,
N.
Ananth
,
F. R.
Manby
, and
T. F.
Miller
 III
, “
Exact nonadditive kinetic potentials for embedded density functional theory
,”
J. Chem. Phys.
133
(
8
),
084103
(
2010
).
21.
D. G.
Jason
,
T. A.
Barnes
, and
T. F.
Miller
 III
, “
Embedded density functional theory for covalently bonded and strongly interacting subsystems
,”
J. Chem. Phys.
134
(
16
),
164108
(
2011
).
22.
S.
Fux
,
C. R.
Jacob
,
J.
Neugebauer
,
L.
Visscher
, and
M.
Reiher
, “
Accurate frozen-density embedding potentials as a first step towards a subsystem description of covalent bonds
,”
J. Chem. Phys.
132
(
16
),
164101
(
2010
).
23.
D.
Schnieders
and
J.
Neugebauer
, “
Accurate embedding through potential reconstruction: A comparison of different strategies
,”
J. Chem. Phys.
149
(
5
),
054103
(
2018
).
24.
Y. G.
Khait
and
M. R.
Hoffmann
, “
On the orthogonality of orbitals in subsystem Kohn–Sham density functional theory
,” in
Annual Reports in Computational Chemistry
(
Elsevier
,
2012
), Vol. 8, pp.
53
70
.
25.
P. K.
Tamukong
,
Y. G.
Khait
, and
M. R.
Hoffmann
, “
Accurate dissociation of chemical bonds using DFT-in-DFT embedding theory with external orbital orthogonality
,”
J. Phys. Chem. A
121
(
1
),
256
264
(
2016
).
26.
H. J. C.
Berendsen
,
D.
van der Spoel
, and
R.
van Drunen
, “
GROMACS: A message-passing parallel molecular dynamics implementation
,”
Comput. Phys. Commun.
91
(
1-3
),
43
56
(
1995
).
27.
J.
Wang
,
R. M.
Wolf
,
J. W.
Caldwell
,
P. A.
Kollman
, and
D. A.
Case
, “
Development and testing of a general amber force field
,”
J. Comput. Chem.
25
(
9
),
1157
1174
(
2004
).
28.
W. L.
Jorgensen
and
J. D.
Madura
, “
Temperature and size dependence for Monte Carlo simulations of TIP4P water
,”
Mol. Phys.
56
(
6
),
1381
1392
(
1985
).
29.
J.
VandeVondele
,
M.
Krack
,
F.
Mohamed
,
M.
Parrinello
,
T.
Chassaing
, and
J.
Hutter
, “
Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach
,”
Comput. Phys. Commun.
167
(
2
),
103
128
(
2005
).
30.
We specifically used the DZVP-GTH-PADE basis for all atoms, in the CONFINED variant where available. The basis-set is available from https://github.com/SINGROUP/pycp2k/blob/master/examples/BASIS_SET.
31.
P.
Giannozzi
,
O.
Andreussi
,
T.
Brumme
,
O.
Bunau
,
M.
Buongiorno Nardelli
,
M.
Calandra
,
R.
Car
,
C.
Cavazzoni
,
D.
Ceresoli
,
M.
Cococcioni
 et al, “
Advanced capabilities for materials modelling with quantum espresso
,”
J. Phys.: Condens. Matter
29
(
46
),
465901
(
2017
).
32.
Y.
Gao
,
D.
Neuhauser
,
R.
Baer
, and
E.
Rabani
, “
Sublinear scaling for time-dependent stochastic density functional theory
,”
J. Chem. Phys.
142
(
3
),
034106
(
2015
).
33.
J.
Towns
,
T.
Cockerill
,
M.
Dahan
,
I.
Foster
,
K.
Gaither
,
A.
Grimshaw
,
V.
Hazlewood
,
S.
Lathrop
,
D.
Lifka
,
G. D.
Peterson
, and
Others
, “
XSEDE: Accelerating scientific discovery
,”
Comput. Sci. Eng.
16
(
5
),
62
74
(
2014
).