Communication: Excited states, dynamic correlation functions and spectral properties from full configuration interaction quantum Monte Carlo

In this communication, we propose a method for obtaining isolated excited states within the full configuration interaction quantum Monte Carlo framework. This method allows for stable sampling with respect to collapse to lower energy states and requires no uncontrolled approximations. In contrast with most previous methods to extract excited state information from quantum Monte Carlo methods, this results from a modification to the underlying propagator, and does not require explicit orthogonalization, analytic continuation, transient estimators, or restriction of the Hilbert space via a trial wavefunction. Furthermore, we show that the propagator can directly yield frequency-domain correlation functions and spectral functions such as the density of states which are difficult to obtain within a traditional quantum Monte Carlo framework. We demonstrate this approach with pilot applications to the neon atom and beryllium dimer.

Almost all of the many variants of projector quantum Monte Carlo (QMC) rely on the properties of the operator e^−βH, where due to its similarity to the time-evolution operator, the variable β is denoted “imaginary time.” Generally, this imaginary time is discretized, and the operator iteratively applied as a short-time propagator, in order to simulate its action in the large β limit.¹ Expressing an initial wavefunction in the eigenbasis of the Hamiltonian of interest, the application of this e^−βH propagator results in a projection onto the

$i th$ eigenstate proportional to

$e^{-\beta E_i}$

$e−βEi$⁠, where E_i is the energy of this eigenstate. It is clear to see that in this large β limit, the projection onto the lowest energy eigenvector dominates the wavefunction, whereby ground state properties can be extracted, assuming some overlap with the initial wavefunction. While this formalism is clearly powerful, by construction it exponentially quickly projects out excited states of the system which may be of interest, and are of critical importance in the simulation of finite-temperature properties, reaction dynamics, photochemistry and many other areas.

To date, isolating excited states of systems via projector QMC methods has only been practical with a restriction on the projection to sample a space which is approximately orthogonal to those of the lower energy states, via nodal constraint,² or orthogonalization against them in a subspace projection method.^3,4 More indirectly, statistical methods have been used on short periods of imaginary-time in order to isolate individual decay rates in the spectrum by analytic continuation to a real-time dynamic.^5,6 However, these approaches are not entirely satisfactory; accurate nodal surfaces for excited states can be difficult to obtain, resulting in a larger fixed node error and potentially transient estimators, while the subspace projection method has limited applicability.⁷ In addition, the Bayesian techniques which rely on maximizing entropy to approximate a notoriously unstable inverse Laplace transform, have difficulty achieving quantitative accuracy within noisy data sets.^8,9 Despite this, there are examples of accurate excited states within the nodal constraint,¹⁰ while maximum entropy techniques are particularly prevalent in solid state calculations to obtain the density of states, often in the case of continuous-time QMC as applied to quantum impurity models and dynamical mean-field theory.^11,12

Here, we take a different approach to the calculation of excited states, within the context of full configuration interaction quantum Monte Carlo (FCIQMC).^13–15 This recently introduced method applies the imaginary-time evolution propagator to a stochastic “walker” representation of the wavefunction expressed in the full space of Slater determinants constructed from a single-particle basis of size M. Although this reintroduces a basis set error compared to those methods operating in real space, it confers various advantages which mitigate this. The discrete basis allows for an efficient walker annihilation algorithm, which can exactly overcome the fermion sign problem in the sampling, provided enough walkers are used.¹⁶ The “initiator” approximation was formulated to maintain a high annihilation rate, and control the growth of noise in a systematically improvable fashion.^14,15 This has allowed far larger systems to be treated at an accuracy comparable to that of FCI (exact diagonalization), within small and systematically improvable error bars. In addition, a semi-stochastic adaptation of the algorithm,¹⁷ as well the introduction of a partial nodal constraint¹⁸ and ideas from quantum chemistry¹⁹ hold promise of increased accuracy and efficiency of the method.

In order to project out a targeted excited state rather than the ground state, we propose the use of a projection operator of the form

For sufficiently large β, this Gaussian propagator will result in the dominant eigenstate being the one closest in energy to the chosen value of the diagonal offset S, termed the shift. In the eigenbasis of H, {|Ψ_i⟩, E_i} with Ψ₀ representing the ground state, and starting from an initial wavefunction

$|ψT⟩$ with S = E_k, it can be seen that the long time propagation results in

We note here that a projector of this form was proposed back in 1983 within continuum QMC approaches,^20,21 although due to sign problems, and significantly larger time step errors resulting from the fact that H² is more singular than H, only one-electron systems were reported, and no modern implementation exists in the literature. This issue of the time step highlights another advantage of working in a finite basis, in that the spectrum is bounded both from below and above. This allows for linearization of the short-time propagator,

without becoming unbound and dominated by very high energy states oscillating in time, and without incurring time step errors in the final wavefunction so long as the time step is less than an upper bound given by

$τ≤2E max −E min$⁠. Repeated application of the short-time propagator therefore results in a “power-method” for states on the interior of the spectrum, rather than at the extrema, with similarities to filter diagonalization.^22,23

Propagation with Eq. (1) leads to a theoretical decay of state j from i as

$e−β2((Ei−S)2−(Ej−S)2)$⁠. In contrast with the ground state propagator, this rate of decay depends on S, with it being advantageous to choose S to be as close to the energy of the state of interest as possible. However, even if S is chosen exactly, the projection of the non-dominant states is slower compared to the ground state propagation. In addition, unless S is chosen exactly, the long-time propagation of the dynamic will result in a continued projection onto a decaying function of all states, including the dominant one. For this reason, the factor of A is introduced into the short-time propagator, such that at convergence, its value can be varied in order to maintain a constant L¹ normalization of the dominant wavefunction and a stable number of walkers. This is analogous to the variation of S in the ground state projection.¹³ 

The differential formulation of the exact dynamic governed by this propagator for a given component of the wavefunction, C_i, can be written as

where ε → τ² as A → 1, and the application of H² has been decomposed by a resolution of the identity over the connecting space of determinants j. This formulation is now amenable to stochastic integration with a discrete walker representation of the determinantal wavefunction coefficients C. As with the ground state projection, there is no unique stochastic algorithm for this dynamic, but the one which we found to be efficient involves a double spawning cycle, which requires little additional overhead compared to the ground state algorithm, and no additional memory requirements. Each iteration, the determinants represented by k are run through, and a spawning attempted to determinant j, in the same fashion as the ground state propagation, but in negative time. This results in a spawning probability to a connected determinant j with a stochastically realised signed amplitude of

$τHjkP(k|j)$⁠, where P(k|j) represents the normalised probability to randomly select symmetry-connected determinant j from k. Successfully spawned walkers are subsequently propagated again in the same iteration with a now forwards-time signed amplitude of

$-\frac{\tau H_{ji}}{P(j|i)}$

$−τHjiP(j|i)$⁠. Care must be taken that for determinant k, the diagonal “death” processes from the first application of H are now interpreted as spawning events, which are also subsequently propagated via H_ij. Each iteration, the factor of A is applied initially as a separate enhancement of the local population of each determinant, with the absolute population on each determinant growing with probability AC_k.

In Fig. 1, we present an illustrative example of the algorithm for the helium dimer in a cc-pVDZ basis, small enough such that the full spectrum of eigenstates can be calculated and the convergence of the method analysed. The value of S was fixed at

$∼40 mE h$ higher than the second excited state, but such that it remained the dominant state in the dynamic, and was subsequently found to be correctly projected out over time. This is despite working in a canonical representation, and starting with a single walker on the Hartree–Fock determinant which had an initial overlap with the ground state of close to one. In order to grow the walkers, A was initially fixed at a value of 1.004, and was varied when a target number of walkers was reached, in common with the procedure for the ground state propagation. The convergence of the projected energy, as defined in Ref. 13, is shown in Fig. 2 for the same simulation, and reflects the decay from the sampled wavefunction of the first excited state.

In order to reliably extend to larger molecular systems, it is worth considering how to transfer the salient elements of the initiator approximation into this new propagation. The basis of this approximation is to attempt to propagate walkers corresponding to wavefunction signal exactly, while walkers judged to be potentially noise are propagated with a truncated Hamiltonian which acts only over the instantaneously occupied subspace.^14,15 This is systematically improvable as the instantaneously occupied subspace grows, or the criterion for walkers corresponding to signal becomes more inclusive. The separation between walkers corresponding to signal and potential noise is not unique. However, it seems sensible to retain the tested feature from ground state propagation that newly spawned walkers on previously unoccupied determinants (i) at the end of an iteration, must have come from an initial determinant (k) which is deemed to have a well-established sign, and therefore a population of walkers above

$n add$⁠. Consequently, all walkers from the application of the first Hamiltonian operator are kept, while the information as to whether C_i is above the initiator criterion is passed through to the annihilation stage of the final set of spawned walkers. No walkers are therefore aborted over the resolution of the identity between H² (determinants j in Eq. (5)).

To test this on a larger system, we study an excited state deep in the spectrum of the 10-electron neon atom in a cc-pVDZ basis, with an energy of approximately

$2.5Eh$ above the ground state. Setting S to equal the CISDTQ energy for the corresponding state, and while remaining in a canonical Hartree–Fock basis and starting from a random distribution of walkers throughout the whole space, we achieved a converged energy of −126.2118(4)

${\textrm {E}_{\textrm {h}}}$

$Eh$⁠, compared to the FCI value of −126.21177

${\textrm {E}_{\textrm {h}}}$

$Eh$⁠. This value is

$4.76{\textrm {mE}_{\textrm {h}}}$

$4.76 mE h$ lower than the initial guess provided by CISDTQ. It would be highly advantageous to develop a robust algorithm for varying S dynamically during the run, as is done for the ground state algorithm. This could be used to maximise the rate of growth of walkers or alternatively minimize A, both of which should adjust S to more closely match the eigenvalue of the state, remove reliance on the initial guess and increase the convergence rate. However, since this requires finding a minimum in a quadratically varying and noisy dataset, no robust algorithm has been identified so far.

Dynamic correlation and response functions due to some perturbation, either in the frequency or time domain, are of critical importance in electronic structure theory,²⁴ and are directly measured in experiments to probe the electronic properties of materials through techniques such as neutron scattering or photoelectron spectroscopy. Many methods, including in general QMC approaches, have significant difficulty in calculating these quantities,²⁵ often having to rely on unstable analytic continuation from imaginary time correlation functions,^5,6,8,9 while other methods can bias towards high or low energy regimes.¹² Although other correlation functions are possible, here we look at the example of an advanced Green's function, a central concept in electronic structure where the “perturbation” at time t = 0 is the creation of a hole in orbital j. For negative time periods, t, these can be written in the time and frequency domain, respectively, as

Unlike the inverse Laplace transform required for the analytic continuation of imaginary time correlation functions, the transform between these two domains is a well-conditioned and numerically stable Fourier transform in the presence of noisy data. Spectral density functions, such as the density of states for extended systems, are then defined in the Lehmann representation²⁶ as

which in the small δ limit tends to

where δ in the above equation represents the dirac-delta function.

Assuming A = 1, application of the propagator in Eq. (1) for a time

$β2=12δ2$ will result in the wavefunction

which when applied to an initial wavefunction

$|ψT⟩=aj|Ψ0⟩$ obtained from the ground-state dynamic, and then projected onto

$\frac{\beta }{\sqrt{\pi }}\langle \Psi _0 |a^{\dagger }_i$

$βπ⟨Ψ0|ai†$ will result in the distribution

for

$S=ω+E0N$⁠. This will tend to the spectral function given in Eq. (10) in the large β limit. The real parts of the Green's function can then be obtained if needed from the Kramers-Kronig relation.²⁴ Results from a pilot investigation of the beryllium dimer in a cc-pVDZ basis, where the exact Green's function can be obtained from complete diagonalization, are shown in Fig. 3.

In order to reduce the statistical error, it may be necessary to average over a small number of independent calculations at each frequency, and this can be combined with an elimination of the bias derived from choosing a correlated sample of

$⟨Ψ0|ai†$ and a_j|Ψ₀⟩,³ by taking the Ψ₀ samples on each side of Eq. (12) from different snapshots in imaginary time. In addition, by storing multiple wavefunctions of the type

$\langle \Psi _0|a^{\dagger }_i$

$⟨Ψ0|ai†$ at the same time, all M² single-particle Green's functions can be calculated at a cost of

$\mathcal {O}[M]$

$O[M]$ FCIQMC calculations per frequency point, without the expectation of any variation in accuracy between high and low energy regimes.

However, despite modest successes, it is clear that obtaining converged results through the use of this operator is substantially more difficult than with the ground state projection. This is mainly due to an additional factor of (τΔE)⁻¹ in the number of iterations required to project out undesired states with energy gap ΔE for comparable accuracy to the ground state propagation. The result is that while in the ground state propagation excited states were filtered relatively quickly with only isolated convergence issues in the case of near degeneracy,¹⁵ the number of iterations required for excited state propagation are substantially increased, as well as the dynamic being less well-conditioned with respect to walker fluctuations. This is also exacerbated by a generally more multiconfigurational wavefunction which increases random error in the projected energy estimator.¹³ A more judicious choice of orbital basis and initial conditions optimized for the state of interest, as well as a multireference projected energy formulation¹⁷ would ameliorate many of these issues. In addition, there is the possibility of preconditioning techniques familiar from iterative diagonalization methods²⁷ being transferred into the stochastic dynamic, as well other operators, such as e^−β|H|, which should behave more efficiently and allow for extension to larger systems. Research in these directions is currently under way. It is clear that alternative propagators within the FCIQMC dynamic holds promise for obtaining accurate excited states.

We would like to acknowledge financial support from the U.S. Department of Energy, Office of Science, through Award Nos. DE-SC0008624 and DE-FG02-07ER46432. We thank Cyrus Umrigar for helpful comments.