Fully quantum embedding with density functional theory for full configuration interaction quantum Monte Carlo

TL;DR

This paper introduces a fully quantum embedding approach for i-FCIQMC using a Huzinaga projection, enabling accurate calculations of large systems like molecules with bonds to benzene, overcoming previous convergence issues.

Contribution

The work develops a novel quantum embedding method for i-FCIQMC that improves convergence and accuracy for large, complex systems, extending its applicability.

Findings

Embedding enables i-FCIQMC to converge to CCSD(T) benchmarks.

Without embedding, i-FCIQMC struggles with large systems.

Embedding reduces computational difficulty and improves accuracy.

Abstract

In common with many high-accuracy electronic structure methods, the initiator adaptation of full configuration interaction quantum Monte Carlo (i-FCIQMC) has difficulty treating realistic systems with large numbers of electrons. This barrier has prevented the application of i-FCIQMC to questions of catalysis that, even for the simplest of models, require high-accuracy modeling of several features of the electronic structure, such as strong and dynamic correlation, and localized vs. delocalized bonding. We here present a fully-quantum embedded version of i-FCIQMC , which we apply to calculate the bond dissociation energy of an ionic bond (LiH) and a covalent bond (HF) physisorbed to a benzene molecule. The embedding is performed using a recently-developed Huzinaga projection operator approach, which affords good synergy with i-FCIQMC by minimizing the number of orbitals in the calculation. We find that, without embedding, i-FCIQMC struggles to converge these calculations due to their substantial system sizes and a lack of error cancellation between reactants and products. With embedding, the i-FCIQMC calculation converges straightforwardly to CCSD(T) benchmarks. Our results suggest that embedded i-FCIQMC will be able treat system sizes well beyond our current reach (even though embedding introduces an error). We discuss how embedding might be improved (and thus the introduced error reduced) using i-FCIQMC energies as benchmarks.

Equations

$$\frac{d|{\Psi }_0⟩}{d\tau }=-\hat{H}|{\Psi }_0⟩$$

$$\underset{\tau \to \infty }{\lim }{e}^{-(\tau \hat{H}-S)}|D_0⟩\propto |{\Psi }_0⟩$$

$$|{\Psi }_0⟩=\sum _i{c}_i|D_i⟩$$

$$c_i^{m+1}-c_i^m=c_i^m\tau (-H_{ii}+S)-\sum _{j\ne i}{c}_j^m\tau H_{ij}.$$

$$H_{ij}=⟨D_i|\hat{H}|D_j⟩.$$

$${\gamma }^{\text{tot}}={\gamma }^{\text{A}}+{\gamma }^{\text{B}}$$

$${\mathbf{F}}^{\text{A-in-B}}={\mathbf{h}}^{\text{A-in-B}}[{\gamma }^{\text{A}},{\gamma }^{\text{B}}]+\mathbf{g}[{\gamma }^{\text{A}}]$$

$${\mathbf{h}}^{\text{A-in-B}}[{\gamma }^{\text{A}},{\gamma }^{\text{B}}]=\mathbf{h}+\mathbf{g}[{\gamma }^{\text{A}}+{\gamma }^{\text{B}}]-\mathbf{g}[{\gamma }^{\text{A}}]+{\mathbf{P}}^{\text{B}}$$

$${\mathbf{P}}^{\text{B}}=-\frac{1}{2}\left({\mathbf{F}}^{\text{AB}}{\gamma }^{\text{B}}{\text{S}}^{\text{BA}}+{\text{S}}^{\text{AB}}{\gamma }^{\text{B}}{\mathbf{F}}^{\text{BA}}\right),$$

$$E_{\mathrm{total}}=E_{\mathrm{KS}-\mathrm{DFT}}^{\mathrm{total}}-E_{\mathrm{DFT}-\mathrm{in}-\mathrm{DFT}}^{\text{A}}+E_{\mathrm{iFCIQMC}-\mathrm{in}-\mathrm{DFT}}^{\text{A}},$$

$$\begin{split}\ce{LiH&->Li++H-}\\ \ce{C6H6\bond{-}LiH&->C6H6\bond{-}Li++H-}\\ \ce{HF&->H++F-}\\ \ce{C6H6\bond{-}HF&->C6H6\bond{-}F-+H+}\end{split}$$

$$E_{\mathrm{corr},\mathrm{A}-\mathrm{B}}=E_{corr,A}-E_{corr,B}+A_1{e}^{-N_w^{{\alpha }_1}}-A_2{e}^{-N_w^{{\alpha }_2}}$$

$$⟨\Psi (\tau )|{\Psi }_i⟩=C_i\exp (-\tau (E_i-E_0))$$

$$\text{ce}C6H6\text{bond}-HF->C6H6\text{bond}-H\text{bond}...F$$

Full text

\altaffiliation

These authors contributed equally to this paper First University] Department of Chemistry, University of Iowa \alsoaffiliation[Another] University of Iowa Informatics Initiative, University of Iowa

\altaffiliationThese authors contributed equally to this paper Second University] Department of Chemistry, University of Minnesota

First University] Department of Chemistry, University of Iowa \alsoaffiliation[Another] University of Iowa Informatics Initiative, University of Iowa

Second University] Department of Chemistry, University of Minnesota

First University] Department of Chemistry, University of Iowa \alsoaffiliation[Another] University of Iowa Informatics Initiative, University of Iowa

1 Introduction

Catalysis often involves bond rearrangements at surfaces, a process featuring closely-separated energy minima, stretched bonds, and transition states. The electronic structure of these systems can become extremely complex; combined with energy differences that can be sub-millihartree, systematic study of catalytic bond rearrangements necessitates the development of new high-accuracy quantum chemistry methods. Although this is a subject of active and ongoing investigation, the high cost of wavefunction methods in particular prevents their widespread application.

One such method is full configuration interaction quantum Monte Carlo (FCIQMC) and its initiator adaptation ($i-$FCIQMC), which are both members of a family of particularly attractive high-accuracy electronic structure methods that seek to combine the exactness of full configuration interaction (FCI) with the increased speed achieved by quantum Monte Carlo (QMC).1, 2 The first FCIQMC paper showed that the FCI ground-state wavefunction could be stochastically sampled due to the sparsity in the Hamiltonian; 1 it had previously been considered that there was no way to sample such a large vector as the exact FCI wavefunction. Since this pioneering work, many adaptions to FCIQMC and $i-$FCIQMC have been developed successfully for calculating correlation energies of a wide variety of benchmark systems.

$i-$FCIQMC has already been used for a variety of applications on relatively small systems, including model systems (such as the Hubbard model3, 4 and the uniform electron gas5, 6) and dimers (such as \ceC27 and \ceCr2 8). It has also seen real applications that are more ambitious, such as iron porphyrins, which used a complete active space adaptation,9 and fully periodic nickel oxide chains 10. A significant amount of investigation has also been aimed at using the full scheme and initiator adaption of FCIQMC to stochastically sample reduced density matrices within the FCIQMC method. 11, 12, 13, 14, 15, 16 Altogether, FCIQMC and its adaptations seem well-poised for answering important questions about the electronic structures of complex chemical systems with high accuracy.

Unfortunately, like all of its high-accuracy cousin methods, $i-$FCIQMC is limited in its scope by its high cost: it can only treat relatively small system sizes (which here means number of electrons). Many further adaptations to $i-$FCIQMC have been developed to allow for the application of $i-$FCIQMC to larger systems. These adaptations include a combination of complete active space self-consistent field (CASSCF) with $i-$FCIQMC 9, the semi-stochastic projector Monte Carlo method 17, model space QMC,18 19 heat-bath configuration interaction,20 perturbation theory 21, stochastic multi-configurational self-consistent field theory (MCSCF) utilizing the FCIQMC methodology 14, use of a transcorrelated Hamiltonian with $i-$FCIQMC 22, 23, and combinations of the above methods, such as semistochastic heat-bath CI 24, 25, 26, 27.

These efforts are made all the more relevant because there are also varieties of FCIQMC which broaden its applicability. Density matrix QMC,28 has been developed for temperature-dependent electronic structure. Additionally, several FCIQMC methods have been developed for use on excited states, such as changing the underlying propagator 29, the Krylov-projected QMC method 30, utilizing a Löwdin partitioning technique18, using a Gram–Schmidt procedure 31, and by restricting the population to the orthogonal complement of the low lying states 32. The stochastic approach in the Slater determinant space has also been studied on the coupled cluster equations, called coupled cluster Monte Carlo. 33, 34, 35, 36 FCIQMC has also been adapted to treat the Clock Hamiltonian, to simulate the full time evolution of a quantum system. 37 A deterministic version of FCIQMC has been developed, 38 as well as a fast randomized iteration framework to essentially perform FCIQMC without walkers. 39 We also note that there are a number of methods which fall under the umbrella of selected configuration interaction (CI), where the CI is solved deterministically, which form a distinct and related family of methods.40, 41, 42

Quantum embedding methods were specifically developed to reduce the problem of scaling present in high level-methods such as $i-$FCIQMC. Embedding methods limit high-level calculations to a small subsystem that is embedded in the potential arising from the rest of the system, reducing the overall computational cost. When highly accurate embedding potentials are used, good accuracy can be achieved even when a subsystem is limited to a few atoms; therefore, embedding methodologies have been successfully applied to a wide variety of systems. 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58 Additionally, a large amount of work has been performed developing accurate embedding frameworks including quantum mechanics / molecular mechanics (QM/MM),59 ONIOM,60 density matrix embedding theory (DMET),61 Green’s function embedding,62, 63 and density functional theory (DFT) embedding.64, 65, 66, 67, 68, 69, 70 A recent review has considered the comparisons between DMET, Green’s function embedding, and DFT embedding and we direct the interested reader to Ref. 71. Many wavefunction methods such as density-matrix renormalization group (DMRG),72 coupled cluster singles and doubles with perturbative triples (CCSD(T)),73 second order Møller-Plesset perturbation theory (MP2),51 and multireference configuration interaction (MRCI)74 have been embedded as the high-level theory; this work presents the first use of $i-$FCIQMC embedding.

The quantum embedding for this work was done using projection-based embedding,73 which is DFT embedding method. Projection-based embedding is one solution to the non-additive kinetic energy problem of DFT embedding.75, 76 The initial projection operator applied to this problem was the $\mu$ projection operator developed by the Manby and Miller groups.73 This projection operator allows two embedded DFT subsystems (DFT-in-DFT) to exactly recreate full-system Kohn-Sham DFT. However, when embedding a wavefunction (WF) subsystem within a DFT environment (WF-in-DFT), the number of orbitals in the WF subsystem is the same as the number of orbitals in the full system. Since WF methods scale poorly with number of orbitals, basis set truncation methods were developed to reduce computational cost.77, 78 The more recent truncation method removes basis functions from a subsystem when the density of that subsystem is below a threshold—a manner that maintains a high degree of accuracy. By decoupling the WF calculation from the total size of the system, WF-level energies may be calculated for systems consisting of hundreds of atoms. The $\mu$ operator method has shown a high degree of accuracy for transition-metal and enzyme catalysis, and oxidation potentials of molecules in solution, among other systems of interest.76, 79 Additionally, several groups have used the $\mu$ projection operator to embed multireference wavefunction methods for application to transition metal catalysts. 80, 81 These systems are inherently multireference; however, as the multireference character is localized to the metal center, $\mu$ embedding calculations were able to closely match experimental results.

Kállay and co-workers introduced the Huzinaga projection operator for DFT embedding;82 however, that work truncated the orbitals by using local correlation methods. We showed that the Huzinaga projection operator could be used for aggressive truncation of the orbital space, where the densities could be absolutely localized on the atomic basis functions centered on atoms within the subsystem.83 This allows for high computational efficiency as the WF subsystem has a greatly reduced number of molecular orbitals. Huzinaga projection embedding has also been successfully extended to periodic systems,84 allowing for cluster or periodic WF calculations embedded in a periodic DFT environment. Given that the absolutely localized basis used in Huzinaga projection-based embedding reduces the number of orbitals to only those centered on the atoms of interest, we here determine the effectiveness of $i-$FCIQMC on a absolutely localized subsystem within the embedding potential of the full system.

We are generally motivated to increase the range and scope of systems available for study with $i-$FCIQMC. With a view toward our long-term interests in the study of bond-breaking and bond rearrangement on surfaces relevant to heterogeneous catalysis, we here study bond dissociation for diatomic molecules containing ionic or covalent bonds (specifically, LiH and HF, respectively) physisorbed onto a benzene molecule using $i-$FCIQMC. This type of calculation (with $\sim 35$ active electrons) is currently at the edge of applicability for $i-$FCIQMC; sometimes the system can be treated, and other times it cannot be treated. We show that embedding greatly alleviates the cost scaling of our model system. Specifically, data show that $i-$FCIQMC calculations performed on the full system (including both the diatomic molecule and the benzene molecule) fails to converge, whereas the system in which the benzene is represented by embedding converges with the same efficiency as an isolated molecular calculation. We analyze the type of convergence behaviors in $i-$FCIQMC and relate them to the differing electronic structures of the dissociation reactants and products. In addition, we explore the applicability of $i-$FCIQMC to a range of atomic separations of HF on benzene by calculating a dissociation curve using both $i-$FCIQMC and CCSD(T). We show that for HF on benzene, $i-$FCIQMC does not have the same failure CCSD(T) shows in regions of strong correlation.

2 Methods

2.1 $i-$FCIQMC

Full configuration interaction quantum Monte Carlo 1 and its initiator adaptation 2 attempt to solve for the ground-state wavefunction $|\Psi_{0}\rangle$ of the imaginary-time Schrödinger equation of a given Hamiltonian $\hat{H}$:

[TABLE]

where $\tau$ represents imaginary time. Beginning with a wavefunction that has non-zero overlap with the ground state, this equation can be solved in the long-imaginary-time limit to give the ground state wavefunction:

[TABLE]

where $|D_{0}\rangle$ is the reference Slater determinant, here taken to be the Hartree–Fock wavefunction. This relationship holds for any constant energy shift $S$. When long enough imaginary time $\tau$ has passed, $S$ can be averaged, and the correlation energy ($E_{\mathrm{corr}}=E_{\mathrm{total}}-E_{\mathrm{Hartree-Fock}}$) found.

The full configuration interaction wavefunction is typically written as a sum of Slater determinants, $|D_{i}\rangle$,

[TABLE]

As such, the imaginary time evolution operator acts in a determinant space.

Substituting Eq. (3) into Eq. (1) gives an expression which can be written as a finite difference

[TABLE]

Here, $c_{i}^{m}$ is the coefficient of the $i^{\mathrm{th}}$ determinant at the $m^{\mathrm{th}}$ iteration of the simulation (after which $m\tau$ units of imaginary time have elapsed). The Hamiltonian is represented in the Slater determinant basis as:

[TABLE]

In the original FCIQMC algorithm, the weight $c_{i}$ takes integer values. 1 The walker population $N_{w}$ is given by $N_{w}=\sum_{i}c_{i}$. When $S$ is varied to keep the walker population constant, its average becomes an estimate of the total ground-state energy.

The population of particles evolves towards the ground state using the following three steps introduced by Booth et al:

The particles with weight $c_{i}$ are allowed to spawn from site $i$ to a connected site $j$, where $H_{ij}\neq 0$ and $i\neq j$. The probability of spawning, $p(j|i)$ is uniform over the $j$ which are connected by one or two electron excitations to $i$. The integer part of $\frac{H_{ij}\tau}{p(j|i)}$ (including its sign) is then added to the weight at $j$. The non-integer remainder $r$ is added with probability $|r|$ as $\pm 1$, where the sign comes from the sign of $r$. 2. 2.

Each particle with weight $c_{i}$ changes its weight by $|S-H_{ii}|\tau$. As above, the integer part of $|S-H_{ii}|\tau$ is added to the weight at $i$. The non-integer remainder $r$ treated as above. 3. 3.

Pairs of particles on the same site with opposite weight $c_{i}$ annihilate each other and and are removed from the simulation, leaving a population containing only a single sign on each site.

FCIQMC is not restricted to using only integer weights $c_{i}$. Real weights can be used; this adds a step to the above algorithm where the real weight is rounded off stochastically below a certain threshhold (here, 0.01), chosen to reduce stochastic error and raise efficiency 17.

The initiator adaption to FCIQMC, $i-$FCIQMC , separates the Slater determinant space into those with $n_{add}$ (here, 3) or more walkers and those with fewer. If the origin of a spawning event (item 1. in the list above above) is not an “initiator” and the spawning is attempted onto a site without walkers, $H_{ij}$ is zeroed. The result is a dynamically-modified Hamiltonian, which profoundly influences convergence of the simulation. A simulation is only converged in the limit when changing the walker population no longer changes the energy (i.e., $N_{w}\rightarrow\infty$). This is an important practical limitation that must be contended with when running an $i-$FCIQMCcalculation, and this ensures the wavefunction must be sampled with sufficient detail in order to attain statistical and systematic convergence. As $N_{w}\rightarrow\infty$, the full configuration interaction (i.e. exact) limit is achieved; away from this limit, the calculation contains a small error termed the initiator error. This error typically converges as $\sim\exp(-\alpha N_{w})$ and is challenging to extrapolate away. Reducing this error is crucial to the success of $i-$FCIQMC ; its pre-factor/rate of decay is highly system dependent, and for larger systems can bottleneck the calculations.

2.2 Embedding

To perform $i-$FCIQMC-in-DFT embedding, the full system density is first split into two subsystems, subsystem A and subsystem B

[TABLE]

where $\gamma^{\text{A}}$ and $\gamma^{\text{B}}$ are the densities matrices of subsystems A and B, respectively. We then obtain the DFT densities of the subsystems through a freeze-and-thaw algorithm.83 This algorithm works by iteratively relaxing the density of subsystem A within the embedding potential and projection operator generated by the frozen density of subsystem B, and then freezing the subsystem A density and relaxing the subsystem B density within the embedding potential and projection operator generated by subsystem A until both subsystem densities have converged. The Fock matrix of subsystem A embedded in subsystem B can be written as

[TABLE]

where $\mathbf{g}$ contains the Coulomb and exchange-correlation potential for DFT—and the embedded core Hamiltonian is

[TABLE]

where $\mathbf{h}$ is the one electron Hamiltonian, and thus contains the kinetic and nuclear potential operators for both subsystems, and $\mathbf{P}^{\text{B}}$ is the Huzinaga projection operator for subsystem A, given by

[TABLE]

where $\mathbf{F}^{\text{AB}}$ and $\mathbf{S}^{\text{AB}}$ are elements of the total Fock matrix and overlap matrix described over the basis functions of subsystems A and B. These equations are then analogously defined for the Fock matrix of B in A. Upon freeze-and-thaw convergence at the DFT level, the $\mathbf{h}^{\text{A-in-B}}[\gamma^{\text{A}},\gamma^{\text{B}}]$ is used as the one-electron Hamiltonian for the $i-$FCIQMC calculation; thus, embedding only influences the one-electron integrals for the $i-$FCIQMC calculation. The final embedding energy is then

[TABLE]

where $E_{\mathrm{KS-DFT}}^{\mathrm{total}}$ is the full-system Kohn-Sham (KS)-DFT energy, $E_{\mathrm{DFT-in-DFT}}^{\text{A}}$ is the DFT energy of subsystem A embedded in the DFT potential of the rest of the system, and $E_{\mathrm{iFCIQMC-in-DFT}}^{\text{A}}$ is the $i-$FCIQMC energy of subsystem A embedded in the DFT potential of the rest of the system.

2.3 Calculation details

The atomic coordinates of the systems under investigation were generated using the dispersion-corrected M06-D3 functional and the aug-cc-pVTZ basis set as implemented in Gaussian16. The geometries are presented in the Supplementary Information. Six frozen orbitals were used for the \ceC6H6\bond-LiH canonical systems, and seven frozen orbitals were used for the \ceC6H6\bond-HF canonical systems. No frozen orbitals were used for the embedded integrals of either system, nor the isolated diatomics. In our implementation, QSoME was modified to output integrals for $i-$FCIQMC using PySCF, 85 which were then read into the HANDE software package.86 For the dissociation curve of \ceHF,MOLPRO 87 was also used taking advantage of an already-existing interface with QSoME. These integrals consisted of single-particle Hartree–Fock eigenvalues ($\epsilon_{i}$) and electron repulsion integrals ($v_{ijkl}$).

The $i-$FCIQMC calculations were performed using the open-source code HANDE-QMC. For the \ceC6H6\bond-LiH system, an imaginary time step of $2\times{10^{-6}}$ a.u. was used with 200,000 reports and 20 Monte Carlo cycles between reports. For the \ceC6H6\bond-HF and \ceC6H6\bond-F^- systems, a smaller time step of $9\times{10^{-7}}$ a.u. was used due to the additional electrons present, with 400,000 reports for the first three target populations and 600,000 reports for the largest three target populations. A larger time step of 0.002 a.u. was used for the isolated LiH, HF and the embedded systems, except for the 5 and 6 Å separations, which used a timestep of 0.0002. In order to converge the calculations with respect to the target population, a range of target populations between ${10^{1}}$ and ${10^{6}}$ was used.

Without the embedding algorithm, the LiH physisorbed on benzene system contains 34 electrons, requires $2.8\times 10^{41}$ determinants, and has a storage cost of 700 MB. After embedding is introduced, the subsystem treated with $i-$FCIQMC is reduced to 4 electrons and $2.9\times 10^{4}$ determinants, with an integral storage cost of 440 KB.

3 Results and discussion

It is common for energy differences to yield better convergence (with respect to excitation rank, for example, in coupled cluster theory) than total energies themselves; this phenomenon, known as error cancellation, is a common benefit of running quantum-chemical calculations. In $i-$FCIQMC (in common with FCIQMC), a walker population of a given size ($N_{w}$) represents the wavefunction. The calculation is only exact if it is converged with respect to this walker number. An under-explored issue of $i-$FCIQMC calculations is that convergence is not faster for energy differences than for individual energies. The dissociation energies of LiH on benzene and HF on benzene represent two paradigmatic examples of how dissociation energies can be extremely challenging and costly to converge in $i-$FCIQMC due to a lack of error cancellation between reactants and products.

We hypothesize that adding benzene to straightforward LiH and HF dissociation energy calculations will cause $i-$FCIQMC to fail in a way that can be remedied by using embedding. To test our hypothesis, we calculate the energy changes associated with four reactions:

[TABLE]

Here, we are required to use the closed-shell ionic dissociation products by the embedding code; an open shell implementation is planned. We note that in \ceC6H6\bond-HF, the H atom is closest to the benzene ring and, following \ceH+ removal and geometry optimization, the \ceF- migrates into the plane of the ring. In particular, we reason that the dissociation energy of a LiH or HF molecule physisorbed to benzene will be significantly more difficult to calculate using $i-$FCIQMC due to non-monotonic energy convergence with system size $N$. In contrast with other methods, $i-$FCIQMC does not show error cancellation between systems that contain different numbers of electrons.

Figure 1 shows data we collected in support of our claim. This data is also presented in table form in the Supplementary Information. Each of these plots is an initiator convergence plot, where the walker population is varied from $10^{1}$ to $10^{6}$, and the energy is computed using $i-$FCIQMC . We plot the $i-$FCIQMC energy differences between reactants and products for the LiH and HF dissociation reactions, and compare these differences to CCSD(T) dissociation energies. CCSD(T) can serve as a good benchmark for initiator convergence: initiator error can vary greatly over many orders of magnitude in energy, and CCSD(T) is generally thought to have systematic error only on the order of 1 millihartree.

Figure 1(a) shows that isolated \ceLiH and \ceHF dissociation energies rapidly converge as a function of walker number, showing complete convergence at $10^{4}$ and $10^{5}$ walkers, respectively. The $i-$FCIQMC and CCSD(T) results are in agreement with each other to within 1 millihartree for $N_{w}\geq 10^{3}$, and within 10 millihartree for the smaller target populations. The \ceHF dissociation converges in an oscillatory manner, because \ceHF is slightly slower to converge than \ceF-; in general, fine-grained oscillatory convergence has been shown in individual calculations. 5 The \ceHF system contains more variability at lower walker numbers than the \ceLiH system, as is expected due to the higher number of electrons present in \ceHF. As we expect, our results show that the isolated systems with small numbers of electrons converge with only modest convergence errors.

In contrast to the isolated molecules, convergence is difficult for the dissociation of molecules physisorbed on benzene. The convergence difficulties for these systems are shown in Fig. 1(b), where the oscillatory behavior observed in Fig. 1(a) is magnified; in the case of HF, we are not able to converge this calculation at all in order to obtain a reaction energy, as the energy difference between $10^{5}$ and $10^{6}$ walkers is approximately -0.0597 hartree. Physisorption onto benzene adds 30 electrons to the isolated molecules; thus, significantly harder convergence is unsurprising. Again, since $i-$FCIQMC does not show error cancellation between systems containing different numbers of electrons, \ceC6H6\bond- HF and \ceC6H6\bond- F- converge at different rates, which causes the energy difference between these two systems to be oscillatory. This is a key result of this manuscript that we explore later in further detail.

In Fig. 1(c), we present the results of the $i-$FCIQMC-in-DFT embedded systems. Since embedding decreases the number of electrons treated directly by $i-$FCIQMC, we are able to converge the $i-$FCIQMC energies of \ceC6H6\bond- LiH and \ceC6H6\bond- HF as easily as isolated LiH and HF. We see similar oscillatory behavior in the embedded calculations as we do for the isolated systems: Target populations $10^{1}$ and $10^{2}$ are still not very accurate. Fortunately, as we increase the target population, we see clear convergence. Comparing the three initiator curves across Fig. 1 reveals a similar convergence trend. This is a very encouraging result, as it shows the $i-$FCIQMC-in-DFT embedding gives convergent results while simultaneously reducing the cost of these calculations significantly.

As computational cost is proportional to walker number, the ability to converge a calculation at $10^{3}$ walkers compared with leaving it unconverged at $10^{6}$ walkers represents a cost savings of at least 1000x. Data we present in the SI additionally show a 1000x savings in memory.

We fully appreciate that there is an unquantified embedding error in these calculations. This causes a change in ordering of the \ceC6H6\bond-HF and \ceC6H6\bond-LiH dissociation energies between Fig. 1(b) and Fig. 1(c) at the CCSD(T) level. For completeness, we note that the difference between CCSD(T) embedded calculations and full-system calculations give us an estimate of the $i-$FCIQMC embedding error as 4.31 millihartree and 8.01 millihartree for LiH and HF, respectively. However, our previous studies have shown that the embedding error can be further decreased by enlarging the wavefunction subsystem.84 Although we are interested in quantifying the $i-$FCIQMC embedding error and using it to benchmark embedded CCSD(T), this analysis is beyond the scope of the proof-of-principle offered by this paper. We now analyze the sources of error and the way that embedding overcomes convergence difficulties in $i-$FCIQMC.

3.1 Analysis of different convergence behaviors in $i-$FCIQMC

There are a number of analyses we can conduct in order to probe the extent of the non-convergent behavior described above in Fig. 1(b)—the case where all electrons in the benzene molecule are fully present in the $i-$FCIQMC calculation. In Fig. 2, the convergence of the reactants and products of dissociation for \ceC6H6\bond-LiH and \ceC6H6\bond-HF are shown. It can be seen from this figure that these calculations are not converged with respect to the number of walkers. This represents a particularly severe case where reactant and product energies actually cross over, which causes the energy differences to oscillate rather than converge smoothly, as observed in Fig. 1(b).

Both \ceC6H6\bond-LiH and \ceC6H6\bond-HF represent different types of challenges in convergence. In \ceC6H6\bond-HF, where reactants and products have the same number of electrons, each $i-$FCIQMC calculation appears to be smoothly converging as a function of walker number. Prior work has established the appearance of such smooth convergence as a stretched exponential in the walker population, $\exp(-N_{w}^{\alpha}),\alpha<<1.0$.7 The decay parameters are highly system-dependent, and as such, two converging calculations could easily cross over one another. The general form of two converging calculations is:

[TABLE]

In the case of HF, the combined initiator error, $A_{1}e^{-N_{w}^{\alpha_{1}}}-A_{2}e^{-N_{w}^{\alpha_{2}}}$, obscures or is much larger than the term $E_{corr,A}-E_{corr,B}$. As a result, the reaction energy fails to converge, instead oscillating even at large walker numbers.

The underlying reason for the differences in convergence between \ceC6H6\bond-HF and \ceC6H6\bond-F- is not known. It seems likely that the form of the stretched exponential is itself related to excited state decays in imaginary time, although this has not been established in the literature. Specifically, the overlap between the simulation wavefunction in imaginary time, $|\Psi(\tau)\rangle$, and the FCI excited states, $|\Psi_{i}\rangle$, is expected to decay exponentially in imaginary time:7

[TABLE]

where $E_{i}$ and $E_{0}$ are the excited state and ground state energy eigenvalues, respectively. In this picture, then, a simulation with insufficient walker population would have to get stuck somewhere between one state and another in a way that cannot be resolved by projecting out over more imaginary time steps, because there is not enough information in each timestep to afford resolution of the ground state.

The case of \ceC6H6\bond-LiH is a little different, since \ceC6H6\bond-Li+ exhibits oscillatory convergence already. This case of oscillatory fine structure has been seen before, such as in studies of the uniform electron gas. 5 This on its own hampers convergence, lending an oscillatory character to the reaction energy independent of whether these calculations are themselves converging to the correct energy.

3.2 Hartree–Fock Population

Another measure by which we can compare the isolated and embedded calculations is the number of walkers present on the Hartree–Fock determinant (shown in Fig. 3). This population is sometimes used as a means to determine convergence of an $i-$FCIQMC calculation, since, in the early phase of an $i-$FCIQMC calculation, it does not vary from its baseline of $[\mathcal{O}(1)]$ walker. The number of walkers on the Hartree–Fock determinant also confirms the different convergence behaviors of the full, isolated, and embedded systems: The embedded and isolated systems have Hartree–Fock populations that grow at the same rate, whereas the full system has many less of these kind of walkers. In terms of the walker dynamics, the larger number of determinants in the full system depletes the signal present on the Hartree–Fock determinant and slows convergence.

3.3 Embedding and the sign problem in $i-$FCIQMC

The sign problem in FCIQMC has been related to the amount of spin frustration in the system: each Slater determinant in the system needs to find its sign over the course of a simulation.4 Specifically, the eigenvalue of a matrix $H_{ij}^{\prime}=\delta_{ij}H_{ij}-(1-\delta_{ij})|H_{ij}|$, where $\delta_{ij}$ is the Kronecker delta, whose eigenstate has entirely non-negative components and contaminates solutions.

The signs in $H$ come from the four-index integrals via the Slater–Condon rules, and so it is important to discuss whether there is a significant change in the integrals due to embedding. In Fig. 4, we show a comparison of two types of integrals that are passed between the embedding code and $i-$FCIQMC for isolated \ceLiH compared with embedded \ceC6H6\bond-LiH. In this case, the LiH eigenvalues are generally lowered by between -0.01 hartree to -0.3 hartree by embedding. The specific ratio for each eigenvalue is plotted against its energy-ordered index in Fig. 4(a), showing that as the eigenvalue becomes higher in energy, it is also affected less by embedding. We show the effect on the electron repulsion integrals in Fig. 4(b), where the distribution of the $\sim 1800$ integrals is presented as a histogram. The molecular orbitals differ between the isolated case and the embedding case and this leads to the small changes in the electron repulsion integrals. From the plots above, we would expect that there is not an increase in the complexity of the sign problem, since most matrix elements remain unchanged.

3.4 Application to bond stretching

Bond dissociation energy curves are frequently used to benchmark new developments in FCIQMC 7. This is in part because CCSD(T) is known to fail due to the strong correlation which occurs as the bond is stretched, leading to certain determinants becoming closer in energy while being strongly coupled.88 In order to highlight the potential benefits of FCIQMC to the study of catalysis we can therefore make comparison between FCIQMC and CCSD(T) for a bond dissociation curve.

Here, we model the dissociation of \ceH-F in \ceC6H6\bond-HF by increasing the \ceH-F bond distance. This represents the following dissociation:

[TABLE]

where the dissociation products, by contrast to Eq. (11), show dissociation by drawing the \ceF- away from the molecule with the rest of the geometry frozen.

Figure 5 shows total and correlation energies calculated at different \ceH-F separations. At above 2Å, the CCSD(T) energy decreases in a manner indicative of strong correlation. By contrast, FCIQMC energies appear to level off to an overall correlation contribution to the bond dissociation energy of approximately $-0.221(1)$ Ha (between the equilibrium separation and 6Å). We note in passing that this is different from the previous correlation contribution to bond dissociation energy of $0.0154(2)$ Ha because these fragments are not the same (see Eq. (14) and Eq. (11)).

4 Conclusions

In summary, we here examined convergence difficulties present when using $i-$FCIQMC to calculate the electronic structure of large systems by exploring the bond dissociations of two prototypical molecules, LiH and HF, physisorbed to benzene. Since $i-$FCIQMC does not show error cancellation between systems with different numbers of electrons, the energy differences between reactants and products tended to oscillate. As a result, dissociation energies calculated from $i-$FCIQMC did not converge. To remedy the convergence issues that $i-$FCIQMC has with large systems, we embedded $i-$FCIQMC in DFT. We showed that this new embedded $i-$FCIQMC was better able to converge dissociation energies, giving results that agree with our CCSD(T) benchmarks. By way of an application, we have also shown the ability of $i-$FCIQMC -in-DFT to more accurately model dissociation curves at atomic separations greater than equilibrium than CCSD(T). This demonstrates the ability of the absolute localization approach for Huzinaga projection-based embedding to treat strongly correlated systems using high-level $i-$FCIQMC wavefunctions embedded in DFT.

Since embedded $i-$FCIQMC also reduces the number of electrons (and thus orbitals) in a calculation, embedded $i-$FCIQMC calculations run with substantially lower cost than full $i-$FCIQMC, alleviating the method’s reduced-exponential cost scaling. Based on our results, we estimate the cost saving to be at least 1000x in compute time and 1000x in memory for the model systems studied here. Whereas for larger systems, $i-$FCIQMC calculations can be computationally intractable while embedded $i-$FCIQMC calculations will remain feasible. There are applications for which CCSD(T) fails to give good answers, such as those involving strong correlation or bond breaking; $i-$FCIQMC can treat these applications with high accuracy. As such, we believe that $i-$FCIQMC emdedded in DFT is a significant and realistic step forward for bringing $i-$FCIQMC towards the routine treatment of real applications, as DFT embedding both alleviates convergence concerns and dramatically reduces the cost of the method.

More broadly, the dissociation curve we calculated represents a situation where strong correlation (bond breaking) was treated by QMC and weak correlation (physisorption) was treated by embedding and DFT; this is likely the best-case scenario for our method. To extend this work, we would move towards real systems. A similar embedding approach as the one we take here has already shown promise for being applied to catalysis.89 For example, quantum embedding was applied to a variety of Co-based catalysts to explore the coupling of the electronic structure of the transition metal to that of the ligand in the hydrogen evolution reaction,90 and to explore the multireference character of these systems that presents challenges for DFT when calculating reaction barriers. 91 Also, Carter and coworkers have used a similar but distinct embedding method to study \ceH2 dissociation on Au and Al nanoparticles.92, 93 In our example of HF dissociation, there was no barrier to dissociation, however in the general case a bond dissociation curve could be used to determine transition state energies and therefore kinetic barrier heights.

We believe that our work is particularly timely because there has been a call from prominent researchers studying oxygen reduction catalysts 94 to focus on the understanding and design of multi-functional active sites for next-generation catalysts and suggested embedding methods could help us get there. In examples such as bifunctional sites or in confinement, we may well expect strong correlation which is quantum mechanically coupled to the environment, and we expect FCIQMC-in-DFT embedding to find applications there.

There are limitations to the embedding approach which are relevant for high-accuracy modelling. The general case where this method will work is when the density from DFT is almost exact. This is especially true for the DFT region, whose density does not change in this method due to our using a frozen density approach. In practice, there are studies that have explored the severity of this approximation95 and found that much can be gained by allowing the errors in the density outside the embedded (here, QMC) region to cancel. It is still very much an open question as to whether a region including strong correlation could be left in the embedded region.

Since strong correlation is often investigated by way of model systems, we note the conditions to apply this approach to model systems as follows. The model system would need to have identifiable localized fragments (here, atoms) and orbitals that are associated with that local fragment which have a dot product rule. The formalism by which the subsystems are divided is exact within a Kohn–Sham formalism and partitions the subsystems to have integer numbers of electrons. In principle it would be possible to use non-integer subsystems which has been applied to 1-D hydrogen chain systems.96, 97

One other limitation of this work is that we have not analyzed the added error in the correlation energy introduced when undertaking embedding, since we believe it is outside of the scope of a proof-of-principle and deserves much more attention on its own. Since CCSD(T) can treat the full systems for the prototypical bond dissociations studied in this manuscript, we could have added a correction to our embedded $i-$FCIQMC arising from the CCSD(T) energy difference between the full and embedded systems; this may be a way forward for future work. It is also of note that the embedding error has been analyzed for CCSD-in-DFT in comparison to CCSD83 and also for CCSD(T)-in-DFT in comparison with experiment;84 we would expect comparable errors at this level of theory. We could also treat a system that is small enough to examine the full system with $i-$FCIQMC , resulting in our being able to benchmark the embedding error for the benefit of other practitioners.

Further work will be forthcoming where this embedding is further developed for excited states and EA/IP calculations; FCIQMC can also be interfaces with other types of calculations such as those with periodic boundary conditions for which a separate periodic code exists. Benchmarking $i-$FCIQMC in comparison to other high-accuracy methods (CCSD(T), DMRG, selected CI) to find the relative advantages and disadvantages of each method represents a very interesting open question. To facilitate this, integral files and output files can be found at https://doi.org/10.25820/data.001111.

In closing, we believe that this study highlights an important step forward for both $i-$FCIQMC and embedding. We believe that the work presented here brings the community one step closer to the routine application of high-accuracy electronic structure to study strongly-correlated systems of chemical and technological interest.

5 Acknowledgements

This research is supported by the Nanoporous Materials Genome Center, funded by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences under Award DE-FG02-17ER16362. DSG and JG acknowledge an award of computer time was provided by the ASCR Leadership Computing Challenge (ALCC) program. This research 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. DSG and JG acknowledge additional computer resources were provided by the Minnesota Supercomputing Institute (MSI) at the University of Minnesota. JJS and HRP acknowledge the University of Iowa for funding and the University of Iowa Informatics Initiative (UI3) for computer resources. The code used throughout this work was HANDE (hande.org.uk), QSoME (github.com/Goodpaster/QSoME), and PySCF (sunqm.github.io/pyscf/).