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.1 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
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,2 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.7 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,10 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.16 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,17 as well the introduction of a partial nodal constraint18 and ideas from quantum chemistry19 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⟩, Ei} with Ψ0 representing the ground state, and starting from an initial wavefunction
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 H2 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
Propagation with Eq. (1) leads to a theoretical decay of state j from i as
The differential formulation of the exact dynamic governed by this propagator for a given component of the wavefunction, Ci, can be written as
where ε → τ2 as A → 1, and the application of H2 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
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
Convergence of the propagation to the second excited state of He2 at 2.5 Å separation in a cc-pVDZ basis. S was fixed at −3.65
Convergence of the propagation to the second excited state of He2 at 2.5 Å separation in a cc-pVDZ basis. S was fixed at −3.65
Convergence of the projected energy estimate to the exact eigenvalue. The reference determinant for the projection was dynamically adjusted to project onto the largest weighted determinant in the sample.
Convergence of the projected energy estimate to the exact eigenvalue. The reference determinant for the projection was dynamically adjusted to project onto the largest weighted determinant in the sample.
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
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
Dynamic correlation and response functions due to some perturbation, either in the frequency or time domain, are of critical importance in electronic structure theory,24 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,25 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.12 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 representation26 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
which when applied to an initial wavefunction
for
High energy window of the spectral function A−(1, 1, ω) as given in Eq. (9) with
High energy window of the spectral function A−(1, 1, ω) as given in Eq. (9) with
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
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)−1 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,15 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.13 A more judicious choice of orbital basis and initial conditions optimized for the state of interest, as well as a multireference projected energy formulation17 would ameliorate many of these issues. In addition, there is the possibility of preconditioning techniques familiar from iterative diagonalization methods27 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.