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.

## I. INTRODUCTION

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)\u2212O(Ne3)$ with the number of electrons *N*_{e}, 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 orbitals^{7} 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.

## II. THEORY

### A. Stochastic DFT

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

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

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

and the one-body Hamiltonian is $H=\u221212\u22072+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\rho $.

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,

where we introduced the number of occupied orbitals (*N*_{occ} = *N*_{e}/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,

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

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

where we abbreviate $\Theta 12\u2261\Theta 12(\mu \u2212H)$ and $\xi \mu \u2261\Theta 12\xi $.

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

where *c*_{i} = ⟨*ϕ*_{i}|*ξ*⟩, so

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

Here, *N*_{s} 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, $\Theta (\mu \u2212H)=12erfc\beta \mu \u2212H$, 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, $\Theta 12=\u2211nan(\mu )Tn(Hscaled)$, where *H*_{scaled} is a scaled Hamiltonian with eigenvalues in the range $\u22121,1$ and *a*_{n}(*μ*) are the Chebyshev coefficients of $12erfc\beta \mu \u2212H12$. Therefore,

where the Chebyshev vectors are obtained recursively, |*ξ*^{n}⟩ = 2*H*_{scaled}|*ξ*^{n−1}⟩ − |*ξ*^{n−2}⟩, and |*ξ*^{n=0}⟩ ≡ |*ξ*⟩.

The Chebyshev expansion makes it possible to analytically determine the chemical potential. Specifically, expand Θ = ∑_{n}*b*_{n}(*μ*)*T*_{n}(*H*_{scaled}), where *b*_{n}(*μ*) are the Chebyshev coefficients of $12erfc\beta \mu \u2212H$. Then, using Eq. (2) gives

where $Rn=Ns\u22121\u2211\xi \u27e8\xi |TnHscaled|\xi \u27e9$. Therefore, *μ* is varied until Eq. (11) is fulfilled.

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

### B. DFT embedding

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

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

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

Therefore,

where $\chi i,\mu =\Theta 12(\mu \u2212H)\chi 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**i*∈

*A*] and then orthogonalize them, to produce $\chi i(r)=\u2211j(S\u221212)ij\varphi j(r)$, where

*S*

_{ij}= ⟨

*ϕ*

_{i}|

*ϕ*

_{j}⟩ is the overlap matrix of the embedded-part atomic orbitals.

The second term in the splitting is associated with *Q* ≡ *I* − *P*, the orthogonal projection to the other (“B”) subspace. Since *Q*^{2} = *Q* and inserting the identity operator $I=|\xi \u2009\u2009\xi |\xi $, we get

where

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,

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, $\Theta \mu \u2212H$, 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, *χ*_{iμ} and $\xi \xaf\mu $, 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 *Qξ* functions, but that orthogonality is lost when we act on *χ*_{i} and on *Qξ* with $\Theta 12(\mu \u2212H)$ in Eqs. (14) and (16). Further note that, as mentioned, the same overall Hamiltonian [and therefore the same $\Theta 12(\mu \u2212H)$] 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 N_{stoch}N_{embed}^{2}, 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.

## III. ALGORITHM

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

Generate

*N*_{s}stochastic orbitals: $\xi (r)=\xb11/d3r$, where*d*^{3}*r*is the volume element associated with the grid. Also create a reasonable initial density*ρ*(), which integrates to the correct number of valence electrons.*r*Determine the one-particle effective potential and Hamiltonian

*H*=*T*+ $v$[*ρ*].For each stochastic orbital, project out the components along the atomic basis functions, i.e., prepare $\xi \xaf=Q\xi $,

where

Determine the correct chemical potential

*μ*as the one that integrates correctly the total density, i.e., from Eq. (11) where now

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

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

Calculate the charge densities for this

*μ*from Eq. (17).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

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 $\Theta 12\mu \u2212H$] 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 $\Theta 12\mu \u2212H$.

## IV. COMPUTATIONAL DETAILS

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.

### A. Structure preparation

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*-nitroaniline^{27} 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.

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.

### B. Stochastic embedding DFT details

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 88^{3} 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\beta \mu \u2212H$ was set at *β* = 0.03 hartree^{−1}, requiring 1173 Chebyshev propagations in acting with $\Theta 12(\mu \u2212H)$.

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

Three sets of calculations were carried out. The first used *N*_{s} = 96 stochastic orbitals, without embedding. The second used the same *N*_{s} = 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 $\Theta 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}

## V. RESULTS

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 *N*_{s} = 96 (second set) is comparable to that with no embedding and *N*_{s} = 256 (third set) and is about 2.5 times the time required to perform calculation with no embedding and *N*_{s} = 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 *N*_{S}. 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 *N*_{S}. In fact, the standard deviation is expected to be proportional to $1/NS$.

. | N_{s} = 96 without
. | N_{s} = 96 with
. | N_{s} = 256 without
. | . |
---|---|---|---|---|

. | embedding . | embedding . | embedding . | Deterministic . |

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 |

. | N_{s} = 96 without
. | N_{s} = 96 with
. | N_{s} = 256 without
. | . |
---|---|---|---|---|

. | embedding . | embedding . | embedding . | Deterministic . |

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.

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.

. | N_{s} = 96 without
. | N_{s} = 96 with
. | N_{s} = 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 |

. | N_{s} = 96 without
. | N_{s} = 96 with
. | N_{s} = 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.

## VI. SUMMARY AND PROSPECTS

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 calculations^{9,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

*H*

_{0}(e.g., that of overlapping fragments). Specifically,

where Θ_{0} ≡ Θ(*μ*_{0} − *H*_{0}) 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.

## ACKNOWLEDGMENTS

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 XSEDE^{33} computational Project No. TG-CHE170058.

### APPENDIX: PROOF THAT EMBEDDING SHOULD REDUCE THE STOCHASTIC ERROR

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} : *i* ∈ *A*} such that it is close to the subspace spanned by states in A,

where

The projected filtered stochastic orbital $\xi \xaf\mu (r)$ is then given by a linear combination,

For states belonging to A, instead of sampling *c*_{i}, we are sampling $(ci\u2212cio)$ with average $(1\u2212cio)\u22480$ and a much smaller standard deviation.