Zombie states are a recently introduced formalism to describe coupled coherent fermionic states that address the fermionic sign problem in a computationally tractable manner. Previously, it has been shown that Zombie states with fractional occupations of spin orbitals obeyed the correct fermionic creation and annihilation algebra and presented results for real-time evolution [D. V. Shalashilin, J. Chem. Phys. 148, 194109 (2018)]. In this work, we extend and build on this formalism by developing efficient algorithms for evaluating the Hamiltonian and other operators between Zombie states and address their normalization. We also show how imaginary time propagation can be used to find the ground state of a system. We also present a biasing method, for setting up a basis set of random Zombie states, that allows much smaller basis sizes to be used while still accurately describing the electronic structure Hamiltonian and its ground state and describe a technique of wave function “cleaning” that removes the contributions of configurations with the wrong number of electrons, improving the accuracy further. We also show how low-lying excited states can be calculated efficiently using a Gram–Schmidt orthogonalization procedure. The proposed algorithm of imaginary time propagation on biased random grids of Zombie states may present an alternative to the existing quantum Monte Carlo methods.

The wavefunction of a multi-electron system is usually accurately described by a linear combination of Slater determinants.1 Each Slater determinant represents a specific configuration—how the available electrons are distributed across spin orbitals. However, in a system with many electrons and many orbitals, the number of possible Slater determinants increases rapidly. Even with sensible approximations and truncation, accurately describing the electronic structure of a system requires a large number of Slater determinants and with that comes computational cost. Introducing an active space of a few most important orbitals and few active electrons is the main method of treating this problem. Active space removes abruptly all configurations outside of it, but even though the contribution of individual configurations outside of the active space can be small, their total contribution can be significant. A smooth transition between active and virtual orbitals in active space methods can potentially be beneficial. Such a transition can be perhaps achieved with quantum Monte Carlo (MC) methods. Diffusion MC2 and Green’s function MC3 propagate walkers in continuum space, which removes the need for a one-electron basis set, but they both fall foul of the fermionic sign problem caused by the antisymmetry of the electronic wavefunction. It has been possible to solve this problem using a fixed node approximation,4 but this has had limited further work. A more robust approach is to introduce a Monte Carlo way of selecting the configurations via a random walk in the manifold of Slater Determinants. Full Configuration Interaction Quantum Monte Carlo (FCIQMC), developed by Alavi and co-workers, works by using a long-time propagation in imaginary time and random walkers to stochastically describe the full configuration interaction wavefunction (FCI).5,6 The method does not converge to the bosonic solution5 and obtains the fermionic ground state without a fixed node approximation. FCIQMC has been successful at obtaining FCI results for large systems that were not previously possible, including the neutral and cationic elements from Li to Mg.5–9 A similar method is Monte Carlo Configuration Interaction (MCCI), developed by Greer,10,11 which originally was applied to the single point energy of water10 and later used to find the dissociation energy of water and Hartree–Fock (HF).12 Like FCIQMC, MCCI uses a Monte Carlo procedure to build a compact wavefunction containing the important configurations intended to give accuracy close to FCI results. Coe et al. continued to develop MCCI and showed that it can generate potential energy curves, approaching levels of chemical accuracy, for a range of small molecules, including N2 and CH4.12 MCCI was then extended to include natural orbitals and second-order perturbation theory, which, particularity, at longer bond lengths was shown to improve the accuracy and convergence time of calculated energies.13 MCCI has also been used to simulate multipole moments achieving good accuracy when compared to FCI values.14 MCCI also showed good accuracy when compared to FCIQMC when calculating ionization energies while only requiring a relatively small number of basis functions when compared to the FCI space.14 The method has also been extended to include state-averaging to allow computation of excited states, which was shown to be effective for both H3 and LiF.15 

However, when used as part of a molecular dynamics simulation, the computational cost of electronic structure calculations, especially Monte Carlo methods, is still significant. Hence, finding methods with acceptable cost of electronic structure calculations, which at the same time would not neglect configurations that normally would be outside of the active space, is of great benefit. In a recent paper,16 we introduced the idea of Zombie States (ZSs), which potentially can combine the ideas of randomness and active space. Zombie states have the possibility of an electron in every orbital in a superposition of “dead” and “alive” states16 so that the occupation of an orbital is fractional. This allows multiple Slater determinant configurations to be described by a single ZS. However, ZSs may not have a well-defined number of electrons such that there is a nonzero probability of the state containing a wrong number of electrons. But, a small superposition of ZSs might be able to describe a wave function with the required number of electrons accurately, and the hope is that a smaller basis set size could be used, reducing the computational cost while maintaining acceptable accuracy.

ZSs can be viewed as fermionic Coherent States (CSs) of a two-level system. Previously, coherent states of the harmonic oscillator have been used to describe bosonic systems in second quantization with the help of Herman-Kluck17 and coupled coherent states propagation methods18 and also with generalized coherent states.19 Second quantization Hamiltonians look like those of coupled oscillators with the difference that the oscillators represent the amplitudes and populations of the orbitals occupied by Bosons. Following the idea of these methods, it seems reasonable to try to treat fermions on a similar level. The complication here is that the elements of Grassmann algebra enter into the standard definition of a fermionic coherent state:20ηj=a(0j+ηj1j), where 1j and 0j are the jth orbital in its occupied and unoccupied vacuum states, respectively; a is a normalization factor; and ηj is an element of Grassmann algebra. The necessity of Grassmann algebra makes computation difficult but is required to maintain the correct permutation antisymmetry of the multi-electron fermionic coherent states as well as the anticommutation of creation and annihilation operators. Zombie states were shown16 to be capable of describing fermionic coherent states while removing the need for Grassmann algebra and the use of Wick’s theorem, which are normally required to evaluate matrix elements between fermionic coherent states.21 

Earlier work16 gave a mathematical treatment of the ZSs’ second quantization for fermions, and creation and annihilation operators were defined with the use of a simple sign change rule, which replaces Grassmann algebra. These results were verified by the reproduction of full configuration interaction1 electronic energies for Li2 and LiH via diagonalization of the electronic structure Hamiltonian in the complete basis of randomly selected ZSs. It was also shown that Zombie states can be used for real time propagation of the electronic wave function. A Fourier transform of real-time Zombie state evolution also has reproduced the exact electronic energy levels in the above-mentioned molecules. Both methods are, however, not very efficient and could be used only as a demonstration that Zombie states’ mathematics was correct. It is therefore necessary to develop better ZS-based methods for finding low-lying state electronic energies in molecules, which are most important in chemistry. In this vein, it is necessary to improve upon the naïve algorithms initially presented to increase efficiency; computation of the Hamiltonian in the original article16 had O(M5) scaling, which can be reduced. This also follows for other key operators. Bosonic coherent states can be defined as an eigenstate of the annihilation operator, whereas the only eigenstate of the annihilation operator for a fermionic coherent state is the vacuum. Thus, we also seek to provide an alternative formalism for Zombie states in the language of second quantization that allows Zombie states to be created from a vacuum state. In this paper, we will present the theoretical advancement in the formalism of Zombie states, deriving a stronger normalization condition. Efficient algorithms for key operators are then given utilizing scaled algorithms and sensible manipulations of the program to greatly reduce the computation time of the Hamiltonian as well as other properties, such as spin and the number of electrons.

Furthermore, imaginary time propagation is validated as an effective method for finding the ground-state energies of states of real fermionic systems, capable of reproducing full CI (in a truncated basis),1 in this case for Li2. We also show how a biasing method can be used to create a small basis of Zombie states, which with imaginary time evolution can find low-lying state energies. To further improve the accuracy, we show how the procedure called “cleaning” removes the contribution of unwanted numbers of electrons and improves the accuracy of electronic energy. We also make an attempt to create a method for finding excited states. It has previously been shown that a Gram–Schmidt procedure can be applied to orthogonalize higher energy states against lower lying ones while adding minimal computational cost.22 We apply a similar Gram–Schmidt procedure to imaginary time propagation of ZSs too, showing that low-lying excited states can be found in this manner.

This article is structured as follows: a brief recap of the original Zombie state formulation is given followed by the algebraic developments; the algorithmic improvements for the two-electron Hamiltonian are detailed; results for imaginary time evolution of Li2 for various basis sets and finally results for imaginary time evolution for excited states of Li2 using Gram–Schmidt orthogonalization are given.

Zombie states16 are coherent antisymmetrized superpositions of “dead” and “alive” electronic states. Considering the single mth spin orbital, a coherent state is

ζ(a1m,a0m)=a1m1m+a0m0m,
(1)

where 1m corresponds to there being an electron in spin orbital m and 0m corresponds to the mth spin orbital being empty, consistent with the conventional electronic structure notation.

We can generalize this to a coherent state, which is a Slater Determinant of one-electron Zombie states,

ζ=ζ1ζ2ζM,
(2)

and which can be summarized by 2M coefficients,

ζ=a11a12a1(m1)a1ma1(m+1)a1Ma01a02a0(m1)a0ma0(m+1)a0M.
(3)

The notation used here is synonymous with that in Ref. 16, where in amjj, mj = 0 refers to the mth spin orbital being “dead” (unoccupied) and mj = 1 refers to the mth spin orbital being “alive” (occupied) and j refers to the spin orbital number. Note that the Zombie state contains 2M complex coefficients, where M is the total number of spin orbitals.

The overlap of two states is

Ωab=ζ(a)ζ(b)=j=1Mmj=0,1amjj(a)*amjj(b).
(4)

Calculation of this is O(M), where there are M spin orbitals, whether or not they are occupied.

The action of the creation and annihilation operators on a ZS is given by

b̂mζ(b)=a11(b)a12(b)a1(m1)(b)a0m(b)a1(m+1)(b)a1M(b)a01(b)a02(b)a0(m1)(b)0a0(m+1)(b)a0M(b),
(5a)
b̂mζ(b)=a11(b)a12(b)a1(m1)(b)0a1(m+1)(b)a1M(b)a01(b)a02(b)a0(m1)(b)a1m(b)a0(m+1)(b)a0M(b).
(5b)

The operators b̂m and b̂m not only act on the mth orbital but also change sign of the amplitudes a1n of alive coefficients for all orbitals with n < m. If a1m and a0n are given by ones and zeros so that Zombie states represent standard Slater determinants, it is easy to see that the sign changing rule is equivalent to the so called Wigner–Jordan rule. As previously shown, all anticommutation relations are as expected.16 Most straightforwardly, but not very efficiently, matrix elements ζ(a)Ĥζ(b) of the second quantized electronic structure Hamiltonian can be computed by sequential application of the creation and annihilation operators,16 

Ĥ=m,nhmnb̂mb̂n+12klmnb̂kb̂lWklnmb̂mb̂n,
(6)

to ζ(b) and overlapping the result with ζ(a),

Hab=ζ(a)Ĥζ(b)=m,nhmnζ(a)ζmn(b)+12klmnWklnmζ(a)ζklmn(b),
(7)

where ζklmn(b)=b̂kb̂lb̂mb̂nζ(b) and ζmn(b)=bmb̂nζmn(b); the overlaps are calculated using Eq. (4).

An electronic wavefunction can be represented as a superposition of K basis Zombie states,

Ψ=kKdkζ(k).
(8)

and the matrix elements described above now allow usage of the ansatz equation (8) for propagation or finding quantum states and their energies. Individual Zombie states generally are not restricted to a particular number of electrons, but for a wavefunction such as Eq. (8), a sufficiently large K number of Zombie states with coefficients dk usually can be chosen such that contributions from unwanted numbers of electrons cancel out. This is very different from standard electronic structure methods, where the problem is to include as many configurations with the right number of electrons as possible.

A standard Hartree–Fock configuration φme(j) that corresponds to me electrons can be written as a Zombie state with “binary” amplitudes of dead and alive states and me ones in the upper row,

φme(j)=100101011010,
(9)

and it can be treated just like any other ZS.

In Subsections II B and II C, algebraic and algorithmic improvements are presented, which allow much more efficient calculation of matrix elements of operators between ZSs and faster propagation of the wave function equation (8). These sections are technical: an alternative derivation of Zombie states from the vacuum state is presented as is the full method we used to reduce Hamiltonian scaling. These are important to the ZS methodology as a whole, but a reader less interested in technical aspects of the method can go directly to the Sec. II D, where it is shown that the right electronic energies can be recovered from imaginary time propagation of the wave function equation (8) with an unrestricted number of electrons.

The previous article16 gave the action of creation and annihilation operators for Zombie states and how the Hamiltonian could be computed for them. Here, we extend the algebraic formalism by showing how a Zombie state can be created from a vacuum state. We show how creation and annihilation operators act in this formalism, as well as computation of the overlap of two Zombie states. Finally, we derive a stronger condition for normalization of a Zombie state.

From the standard electronic structure theory,1 a given electronic occupancy can be written as an antisymmetrized Slater determinant as

φ=m1m2mM,
(10)

where {mi} = 0 if the ith orbital is empty and 1 if occupied, as detailed in Ref. 23 or, equivalently, in the form of Eq. (9). This is equivalent to stating in second quantization notation,

φ=koccb̂k,
(11)

where is the vacuum state, ensuring here and in what follows that the creation operators are applied in the reverse order they appear in the Slater determinant. We could also trivially rewrite Eq. (11) as

φ=j=1M[(1mj)Î+mjb̂j].
(12)

We now consider a generalized form of Eq. (12),

ζ=j=1M(a0jÎ+a1jb̂j),
(13)

where a0j and a1j are complex scalar coefficients. Clearly, if a1j = mj and a0j = (1 − mj), then the same Slater determinant electron state is produced as in Eq. (12).

By inference from Eq. (13), we can define a Zombie operator

ẑja0jÎ+a1jb̂j,
(14)

whose adjoint is

ẑja0j*Î+a1j*b̂j,
(15)

and commutators and anti-commutators

[ẑj,ẑk]=2b̂jb̂ka1ja1k,
(16a)
[ẑj,ẑk]=2b̂jb̂ka1j*a1k*,
(16b)
{ẑj,ẑk}=2(a0ja0kÎ+a1ja0kb̂j+a0ja1kb̂k),
(16c)
{ẑj,ẑk}=2(a0j*a0kÎ+a1j*a0kb̂j+a0j*a1kb̂k)
(16d)

such that Eq. (13) can be written more concisely as

ζ=j=1Mẑj.
(17)

Note that ẑj is a function of two complex numbers, so we can write ẑjẑj(a0j,a1j). This shorthand will become useful later.

The idea of defining new operators from linear combinations of creation and annihilation operators is, of course, nothing new, the most famous example perhaps being Majorana Fermions, which are their own antiparticle, defined as

γ1=12(b̂+b̂),
(18a)
γ2=12i(b̂b̂).
(18b)

However, Zombie operators are not their own antiparticles. We now consider the action of creation and annihilation operators on Zombie operators.

1. Creation and annihilation operators

The action of creation and annihilation operators on a single Zombie operator is given, for jk, by

b̂jẑk(a0k,a1k)=b̂j(a0kÎ+a1kb̂k)=(a0kÎa1kb̂k)b̂j=ẑk(a0k,a1k)b̂j,
(19a)
b̂jẑk(a0k,a1k)=b̂j(a0kÎ+a1kb̂k)=(a0kÎa1kb̂k)b̂j=ẑk(a0k,a1k)b̂j,
(19b)

and if j = k,

b̂jẑj(a0j,a1j)=b̂j(a0jÎ+a1jb̂j)=a0jbja1k(1b̂jb̂j)=ẑj(a0j,a1j)b̂j+ẑj(a1j,0),
(20)
b̂jẑj(a0j,a1j)=b̂j(a0jÎ+a1jb̂j)=a0jb̂j=ẑj(0,a0j),
(21)

where we use the standard relations1 {b̂j,b̂k}=0, {b̂j,b̂k}=0, and {b̂j,b̂k}=δjk.

Using these results, we can compute the action of the creation and annihilation operators on a Zombie state, which is defined by the product of Zombie operators as in Eq. (17),

b̂jζ=b̂jk=1Mẑk(a0k,a1k)
(22a)
=k=1j1ẑk(a0k,a1k)ẑj(a0j,a1j)b̂j+ẑj(a1j,0)×k=j+1Mẑk(a0k,a1k)
(22b)
=k=1j1ẑk(a0k,a1k)ẑj(a1j,0)k=j+1Mẑk(a0k,a1k),
(22c)

where the ẑj(a0j,a1j)b̂j term vanishes because b̂j tries to destroy an electron that is not there.

It also follows that

b̂jζ=b̂jk=1Mẑk(a0k,a1k)
(23a)
=k=1j1ẑk(a0k,a1k)ẑj(0,a0j)k=j+1Mẑk(a0k,a1k).
(23b)

However, Eqs. (22c) and (23b) themselves describe Zombie states with changed coefficients to those of ζ. The Zombie states they refer to agree exactly with the sign-changing rules presented in the original Zombie states paper, though through more convoluted algebra.

These results can equivalently be stated in terms of the {anjj} coefficients. As ζ is a function of the Zombie coefficients {anjj}=a, we can write ζζ(a). Then, let ζ(j)(a(j))=b̂jζ(a). We then see that

a0k(j)=a0k,kj,a1k,k=j,
(24a)
a1k(j)=a1k,k<j,0,k=j,a1k,k>j.
(24b)

If ζ(j)(a(j))=b̂jζ(a),

a0k(j)=a0k,kj,0,k=j,
(25a)
a1k(j)=a1k,k<j,a0k,k=j,a1k,k>j,
(25b)

again in exact agreement with the earlier results.

2. Overlap of two Zombie states

We can use this formalism to compute the overlap of two Zombie states, ζ(a)ζ(b), where both the bra and ket are formed from the action of M Zombie operators as in Eq. (17). We consider a computation recursively, noting that the vacuum state is normalized, namely, =1, and defining

ζj(b)=k=jMẑk
(26)

such that

ζ1(b)ζ(b),
(27a)
ζM+1(b).
(27b)

We, therefore, have

ζj(a)ζj(b)=ζj+1(a)ẑj(a)*ẑj(b)ζj+1(b)
(28a)
=ζj+1(a)a0j(a)*a0j(b)Î+a0j(a)*a1j(b)b̂j+a1j(a)*a0j(b)b̂j+a1j(a)*a1j(b)b̂jb̂jζj+1(b)
(28b)
=(a0j(a)*a0j(b)+a1j(a)*a1j(b))ζj+1(a)ζj+1(b).
(28c)

The second term in Eq. (28b) is of the form ζj+1(a)b̂jζj+1(b), and since neither ζj+1(a) nor ζj+1(b) contain any contributions from electron j, b̂j acting to the left will try to destroy an electron that is not there, giving zero. Similar arguments show why the third term in Eq. (28b) vanishes and why the fourth term survives.

Combining Eqs. (26)(28) gives

ζ(a)ζ(b)=j=1M(a0j(a)*a0j(b)+a1j(a)*a1j(b)),
(29)

which is identical to Eq. (4), which was obtained by different algebraic means in the original Zombie paper.16 

3. Normalization

Clearly, if we wish a Zombie state (the product of M Zombie operators) to be normalized, then the RHS of Eq. (29) must equal unity when ζ(a)=ζ(b). However, we may want to consider a more general case, where we may not necessarily have wavefunctions formed from a fixed M number of Zombie operators but a variable number of Zombie operators.

Algebraically, if we require

j=1M|a0j(a)|2+|a1j(a)|2=1,
(30)

this will ensure that ζ(a)ζ(a)ζ1(a)ζ1(a)=1, but what if we further require

ζj(a)ζj(a)=1
(31)

for all j, 1 ≤ jM.

This problem can also be solved recursively, since trivially ζM+1(a)ζM+1(a)=1. For ζM(a)ζM(a)=1, from Eq. (28), it is required that |a0M(a)|2+|a1M(a)|2=1. However, from Eq. (26), ζj(a)ζj(a)=(|a0j(a)|2+|a1j(a)|2)ζj+1(a)ζj+1(a), and so for Eq. (31) to be true for all j, we require

|a0j(a)|2+|a1j(a)|2=1j.
(32)

Note that Eq. (32) is a stronger criterion than Eq. (30), i.e., satisfying Eq. (32) will ensure that Eq. (30) is satisfied, but not the other way around.

Equation (32) might correspond to normalization of a “one dimensional” ZS. This normalization is equivalent to the normalization used when introducing quantum superposition sampling to the multiconfigurational Ehrenfest method presented in Ref. 24.

In the event that the Zombie state is real (useful for imaginary time propagation to find stationary states), then a simple choice to satisfy Eq. (32) is

a0j(a)=cos(θj),a1j(a)=sin(θj),
(33)

where θj is a real number such that 0 ≤ θj < 2π. This means that a real, normalized Zombie state can be completely described by M real {θj} values. This can be generalized to complex states24 with

a0j(a)=cos(θj),a1j(a)=sin(θj)eiϕj,
(34)

where ϕj is drawn from the interval 0 ≤ ϕj < 2π.

Here, we consider how to efficiently compute the action of the Hamiltonian and other operators on a Zombie state. We show how the electronic Hamiltonian matrix elements with Zombie states ζ(a)Ĥζ(b) can be calculated more efficiently. First, we show how a more efficient algorithm with a smaller prefactor can calculate matrix elements faster, and second, we will present an algorithm with O(M4) scaling as opposed to the O(M5) scaling of the naïve algorithm used in the initial paper16 described above.

1. Reduced prefactor Hamiltonian algorithms

The bottleneck in the computation of FCI states and energies, in general, is evaluation of the two-electron Hamiltonian

Ĥ2=12i,j,k,lijklb̂ib̂jb̂lb̂k.
(35)

It was previously shown16 that the two-electron part of Eq. (6) can be evaluated as

ζ(a)Ĥ2ζ(b)=12ijklMijklζ(a)ζijlk(b),
(36)

where ζijlk(b)=bib̂jb̂lb̂kζ(b). This can be solved with iterative applications of the creation and annihilation operators, meaning that the original algorithm, therefore, requires M4 terms to be evaluated to solve Eq. (36). Each term has four creation/annihilation operators [each of which is O(M) due to the sign changing rule] and one overlap computation, which is O(M) as

Ωab=ζ(a)ζ(b)=m=1Mnm=0,1anmm(a)*anmm(b).
(37)

Consequently, naïve evaluation of Eq. (6) requires M4 creation/annihilation operations and is O(M5) overall.

Many of the two-electron integrals are zero by spin symmetry, i.e., ijkl=0 unless σ(i) = σ(k) and σ(j) = σ(l), where σ(i) is the spin of spin orbital i. One method for accelerating the naïve algorithm is by evaluating only ζijlk(b)=bib̂jb̂lb̂kζ(b) if ijkl0 on spin symmetry grounds. When applied to the restricted Hartree–Fock (RHF) calculation, this simple code modification gave nearly a four times speed up when calculating a matrix element between two ten-orbital ZSs. This is roughly to be expected since three quarters of the two-electron integrals are zero from spin symmetry.

Furthermore, it is also possible to loop only over combinations of {ijkl} already known to be nonzero by altering the loops over orbitals such that only nonzero (by spin) ζijlk(b) are considered in the first place. This produces a comparable, albeit slightly larger time improvement, ignoring zero two-electron integrals.

In an attempt, to accelerate the code further, we note that computing the action of M4 creation and annihilation operators is expensive and look to evaluate fewer of them. Re-examining Eqs. (35) and (36), we see

ζ(a)Ĥ2ζ(b)=12ijklMijklζij(a)ζlk(b),
(38)

where

ζij(a)=ζ(a)b̂ib̂j=b̂jb̂iζ(a),
(39)
ζlk(b)=b̂lb̂kζ(b).
(40)

We can then compute {ζij(a)}i,j and {ζlk(b)}l,k, which requires M2 creation and annihilation operator applications, while overall the Hamiltonian evaluation is O(M5). Having precomputed {ζij(a)} and {ζlk(b)}, we can then evaluate Eq. (38).

From the definition of the overlap integral in Eq. (4), we see that if both the dead and alive Zombie amplitudes are zero for a given state, i.e., for the mth spin orbital,

ζ(c)=a11(c)a12(c)a1(m1)(c)0a1(m+1)(c)a1M(c)a01(b)a02(c)a0(m1)(c)0a1(m+1)(c)a0M(c),
(41)

then the overlap of this state with any other state will be zero. Furthermore, since the Zombie creation and annihilation operators only move the dead and alive amplitudes within a spin orbital, and not between spin orbitals, ζ(c) will continue to have zero overlap with any other state even after application of creation and annihilation operators. Consequently, if a state such as ζ(c) is itself generated by application of creation and annihilation operators, it will not contribute to the two-electron Hamiltonian and so does not need to be considered in further calculations.

We know that16 

b̂ζb̂a1a0=0a1.
(42)

Consequently, if a1 = 0, application of b̂ returns a vanishing state that will have no overlap with any other. This is a rephrasing, in Zombie terms, of the well-known result from electronic structure theory that trying to annihilate an electron, which is not there will return zero. Applied to calculation of the two-electron Hamiltonian, let us consider ζij(a)=b̂jb̂iζ(a). If a1i(a)=0, then b̂iζ(a) returns a vanishing state and will contribute nothing to the two-electron calculation, even after application of b̂j. If b̂iζ(a)ζi(a) is nonzero, then if a1j(a)=0, the Zombie state vanishes. Consequently, when looping over orbital indices in Eq. (38), if the application of a given annihilation operator would return a vanishing state (which can be determined from a1i(a)=0) and similarly for j, l, k, then the all further calculation for that creation/annihilation operator can be discarded.

In addition, considering the overlap calculation equation (4) and rewriting it as

Ωab=ζ(a)ζ(b)=m=1Ma0m(a)*a0m(b)+a1m(a)*a1m(b),
(43)

we see that if any of the a0m(a)*a0m(b)+a1m(a)*a1m(b) terms are zero, then Ωab = 0. This means that calculation can be terminated as soon as one of the terms is zero and the overlap returned as zero, without having to compute any remaining terms. These improvements gave a similar time improvement to previous changes when working with ten orbitals; however, the time improvement was approximately half as good when calculating a matrix element with 50 orbitals. Full details of the algorithmic code racing can be found in Sec. 1.1 of the supplementary material. All of these algorithms still have (approximate) fifth order scaling, i.e., O(M5), although with a much lower prefactor than earlier versions.

2. Lower-scaling Hamiltonian algorithm

It is, however, possible to reduce the scaling of the two-electron Hamiltonian from O(M5) to O(M4). We first note that evaluation of the two-electron Hamiltonian could be written as

ζ(a)Ĥ2ζ(b)=12ijklMijklζij(a)b̂lζk(b).
(44)

We have previously shown that ζij(a) can be calculated in O(M3) and ζk(b) in O(M2), and these can be calculated separately. Using previous methods, even with the calculation improvements, to evaluate ζij(a)b̂lζk(b)i,j,k,l requires O(M3) for looping over i, j, k, then an O(M) for summing over l, and then another O(M) for evaluating the overlap, making O(M5) overall. We, therefore, consider how to evaluate ζij(a)b̂lζk(b) for a specific set of i, j, k, but over all values of l, in O(M) operations, such that the overall algorithm would be O(M4).

We start by simplifying our notation to ζ(a)b̂lζ(b), which yields

ζ(a)b̂lζ(b)=i=1l1a0i(a)*a0i(b)a1i(a)*a1i(b)a0l(a)*a1l(b)j=l+1Ma0j(a)*a0j(b)+a1j(a)*a1j(b).
(45)

Note the minus sign in the first product caused by the action of the annihilation operator. We then define

si=a0i(a)*a1i(b),
(46a)
ei=a0i(a)*a0i(b)a1i(a)*a1i(b),
(46b)
fi=a0i(a)*a0i(b)+a1i(a)*a1i(b).
(46c)

These can all trivially be calculated in O(M) steps. We also define

gl=i=1lei,
(47a)
hl=i=lMfi,
(47b)

which can be calculated recursively in O(M) steps. Then,

ζ(a)b̂lζ(b)=gl1slhl+1,
(48)

which can also be found in O(M) steps.

The reduced scaling improves calculation of two-electron matrix elements over 20 times for ten orbitals calculation and over 100 times with 50 orbital calculations compared to the naïve calculation. Further details are given in Table 4 in the supplementary material. In the case of 50 orbitals, the scaled algorithm does incur a larger peak memory usage, requiring ∼414% more memory for the calculation. However, the fastest algorithm still only needs 2 MB of memory and so the speed increase is worth the relatively larger memory requirement.

We have also developed an algorithm to calculate the number operator, which scales linearly with the number of orbitals in the system. Using a similar method, we also show how to reduce the scaling of the algorithms used to calculate Sz, Ŝz2 and total spin. Details of these two algorithms are presented in  Appendixes A and  B.

Imaginary time propagation, i.e., solution of the Schrödinger equation in imaginary time,

ddβΨ=ĤΨ,
(49)

is commonly used to find the ground state of a quantum mechanical system.5,6,22 Many methods of quantum Monte Carlo exploit it in electronic structure theory. For example, a random walk in the Fock space of Slater determinants has been shown to be a beautiful technique to solve Eq. (49) and to avoid the wave function node problem, which otherwise makes quantum Monte Carlo for fermions difficult. In this section, we will show that a small randomly selected set of Zombie states can provide an accurate description of the ground state of a molecule.

1. Imaginary time algorithm

Consider a wavefunction Ψ that can be expanded in the basis of K Zombie states ζ(k), where K ≤ 2M, and which are normalized but not necessarily orthogonal,

Ψ=kKdkζ(k).
(50)

Then, applying the identity matrix from Ref. 25 for a nonorthogonal basis set l,kζl)Ωlk1ζ(k) and using the reasoning from Ref. 26, we yield

kdkβζ(l)ζ(k)=kζ(l)Ĥζ(k)dk,
(51)

which can be rearranged in matrix notation as

ḋ=Ω1Hd.
(52)

Note how the inverse overlap matrix with the elements Ωlk=ζ(l)ζ(k) is required and how the algebra developed for bosonic coherent states is being applied to fermionic ZSs.26 

2. Cleaning

The propagation of Eq. (52) relaxes the wave function to the lowest energy state of the system.27 As occupations of ZS spin orbitals are fractional Zombie states, their linear combinations are not restricted to a particular number of electrons. However, they can be projected onto a Fock Space with a given number of electrons. The identity operator can be written as a sum of me electron identities Îme,

Î=me=0,MÎme,
(53)

where the identity covering the Fock Space with me electrons is

Îme=jφme(j)φme(j),
(54)

where the sum in j is over all Fock configurations φme(j), with me being number of electrons, which can be viewed as Zombie states with “binary” amplitudes 1 and 0 and me being unit amplitudes of alive electrons on me occupied spin orbitals.

Then, the superposition of Zombie states in Eq. (8) or Eq. (50) can be projected onto Fock Space of me electrons as

Ψme=ÎmeΨ=jcjφme(j),
(55)

where

cj=k=1,Nzsdkφ(j)ζ(k).
(56)

The energy of the me contribution and its norm can then be calculated as

Eme=jici*cjφme(i)Ĥφme(j)
(57)

and

Nme=jcj*cj,
(58)

respectively. The sum of energies of all me contributions is equal to the energy of the whole Zombie wave function Eq. (8) and all me norms add up to 1, i.e.,

ΨĤΨ=meEme
(59)

and

meNme=1.
(60)

For a completely converged Zombie wave function, the energy and norm of all me other than the right number of electrons ne will be zero, but if convergence is incomplete, the energy of the converged me wave function Ψc can be estimated as

ΨcĤΨcEme/Nme
(61)

by division of the me energy by the me norm. Cleaning is effectively the same as the projection-after-variation method presented in Refs. 28 and 29. The only difference being the ZS wave function is multiconfigurational opposed to a single determinant function.

3. Biased bases

In this section, we will consider various ways of how a basis set of ZSs can be chosen to solve Eq. (52). The best choice (in theory) would be a complete set of 2M Slater determinants φme(j) covering all possible electron occupancies from 0 to M electrons in M spin orbitals. Diagonalization of the Hamiltonian matrix in this basis gives all possible energy levels as it contains all possible combinations, from 0 to M electrons, placed in M spin orbitals. A basis of randomly selected 2M ZSs ζ(k) is also complete and is capable of yielding correct quantum energies and wave functions. However, for efficient calculation, the basis set must be as small as possible. We employ a biased basis that sits between the two extremes of complete order and randomness. In complete active space self-consistent field (CASSCF) calculations,30 electrons and orbitals are split into three groups: inactive, active, and virtual. From this, we have designed a biasing method to set up a ZS basis. Inactive, or core electrons, are low lying and are always (at least in a CASSCF calculation) to be occupied; the active orbitals can either be occupied or unoccupied and virtual orbitals are always empty. We take this idea of splitting electrons into different groups and then exploit the fractional occupation of a Zombie state. This allows orbitals to be set as completely “alive” (or “dead”) while also biasing orbitals to be more or less “dead” or “alive.” It should be noted that though we are using the ideas of splitting electrons into different types from CASSCF, this does not make ZS an SCF method any more than using HF energies as a starting point. Amplitudes for Zombie states can be randomly chosen from normal distributions centered around either completely “alive” or completely “dead” values and the distributions’ width set according to its activity level. To satisfy the normalization condition given by Eq. (32), in this paper, only one Zombie amplitude (dead or alive) is generated using the normal distribution and the other is set using Eq. (33). Hence, the low-lying core electrons have “alive” amplitudes set to one, although a very narrow Gaussian could also be employed, which corresponds to a normal distribution width of zero giving a δ function with a1j = 1 and a0j = 0. The active electron normal distributions can be centered to favor either completely alive or dead amplitudes; the point where centering changes from “alive” to “dead” naturally changes with each system. Thinner distributions are used for orbitals at either end of the active electron/orbital group and wider distributions are used for orbitals where the chance, or not, of occupation is relatively equal. Virtual orbitals would have “dead” coefficients set to 1. However, again, narrow distribution near completely “dead” amplitudes can be employed.

4. Excited states

As previously stated, the aim of developing Zombie states is so it can be implemented in no-adiabatic molecular dynamics simulations and so it would be useful to compute low-lying excited states with the low computational expense that Zombie states could provide. It has previously been shown that a Gram–Schmidt orthogonalization procedure can be combined with full configuration interaction quantum Monte Carlo22 and was found to add little extra computational cost while allowing multiple low-energy states to be studied. We will apply the same process here combining imaginary time propagation and Gram–Schmidt. Generally, Gram–Schmidt orthogonalization, for any set of vectors, is defined by

uk=vkj=1k1projuj(vk)
(62)

with the projection operator

proju(v)=u,vu,u.
(63)

In this implementation, the k lowest states of interest are arbitrarily set, creating dk vectors that are not necessarily orthogonal. The vectors are then each orthogonalized,

dk=dkj=1k1projdj(dk),
(64)

where dk denotes an orthogonalized vector, and they are then normalized. Each state dk is then put into the differential equation equivalent to Eq. (52),

ḋk=Ω1Hdk.
(65)

A single time step is taken, and a new set of dk(β+Δβ)=dk(β)+ḋk(β)Δβ is then found. The process then repeats: vectors are orthogonalized using the Gram–Schmidt process and normalized and a time step is taken.

Here, we will demonstrate that imaginary time propagation is an effective and efficient method for finding ground state energies using a Slater determinant basis as a point of reference. Earlier work16 used Li2 to verify the theoretical basis of Zombie states, and so here we too will use Li2 as an example system. PyScf31 (version 1.7.6) using a 6-31G** basis and the Dh symmetry group has been used to calculate one and two-electron integrals. These are then converted from spatial to spin orbitals, in this case, five spatial molecular orbitals (and therefore ten spin orbitals). A Hamiltonian for the system is computed, using the Hamiltonian detailed in Sec. II C 2, with the specified ZSs. The full code used can be found in Ref. 32 All calculations are carried out in atomic units so energies are in Hartrees.

The first trivial case considered is a complete basis of Slater determinant φme(j) Zombie states. Li2 consists of ten spin orbitals, which gives a full basis a size of 210 = 1024 basis functions. In Fig. 1, the six-electron restricted Hartree–Fock (RHF) determinant is used as a starting point and imaginary time evolution gives the neutral ground state energy of Li2. A full sized basis of 210 = 1024 randomly generated basis functions ζ(k) was also tested. Random Zombie states are generated by randomly generating θ, between 0 and 2π, for each orbital, and the dead and alive coefficients are then found by calculating a1 = cos(θ) and a0 = sin(θ), respectively. The random basis is made up of Zombie states consisting of superposition of “dead” and “alive” electrons compared to the Slater determinant basis, which is entirely binary. The initial vector is a superposition of random Zombie states, which is chosen to be equal to the RHF determinant, and that propagation of this using imaginary time evolution (solid line in Fig. 1) matches the result obtained by using the Slater Determinant basis. The final energy of each state obtained through imaginary time propagation is verified by comparison to the eigenvalues found by diagonalizing the complete Slater determinant Hamiltonian.

FIG. 1.

Imaginary time propagation starting with the six-electron restricted Hartree–Fock Slater determinant using a complete Slater determinant basis (dashed red line) and a complete random basis (solid orange line). Both results tend to neutral ground sate energy of Li2 found by diagonalizing the Slater determinant Hamiltonian, shown as a gray dotted line.

FIG. 1.

Imaginary time propagation starting with the six-electron restricted Hartree–Fock Slater determinant using a complete Slater determinant basis (dashed red line) and a complete random basis (solid orange line). Both results tend to neutral ground sate energy of Li2 found by diagonalizing the Slater determinant Hamiltonian, shown as a gray dotted line.

Close modal

When a reduced basis of 200 random Zombie states is used as seen in Fig. 2, the basis set is no longer capable of reproducing the ground state energy. The reduced basis is created by selecting 200 basis functions ζ(k) at random from the complete random basis previously used. The energy obtained via imaginary time propagation, in this smaller basis, no longer matches the value found by diagonalizing the Slater determinant basis. However, the cleaned energy provides a much better estimate of the ground state energy. For a poorly selected basis, the distribution of energy and norm over possible electronic numbers me is shown in Fig. 3, which demonstrates that only a fraction (≈17.5%) of the final wave function belongs to the right Fock space on ne = 6 electrons. The highest fraction of the norm (≈40%) belongs to the 7e Fock space. The cleaned energy E(6e)/N(6e) = −14.598 267 is somewhat above the FCI energy of the 6e ground state (−14.871 914), but the cleaned energy of the 7e Fock Space E(7e)/N(7e) = −14.750 162 is approaching the 7e FCI (−14.858 062) energy. Thus, the randomly selected basis has, by chance, been set up in favor of the anion.

FIG. 2.

Imaginary time propagation starting with the six-electron restricted Hartree–Fock Slater determinant and using a reduced random basis of 200 Zombie sates. The neutral ground sate energy of Li2 found by diagonalizing the Slater determinant Hamiltonian is also shown in gray.

FIG. 2.

Imaginary time propagation starting with the six-electron restricted Hartree–Fock Slater determinant and using a reduced random basis of 200 Zombie sates. The neutral ground sate energy of Li2 found by diagonalizing the Slater determinant Hamiltonian is also shown in gray.

Close modal
FIG. 3.

Energy (top) and norm (bottom) distribution for each number of electrons for a basis set of 200 randomly generated Zombie states. Only a fraction of the final wave function belongs to the Fock space on ne = 6 electrons.

FIG. 3.

Energy (top) and norm (bottom) distribution for each number of electrons for a basis set of 200 randomly generated Zombie states. Only a fraction of the final wave function belongs to the Fock space on ne = 6 electrons.

Close modal

The electronic structure of the Li2 is not completely random and unknown, so it is sensible to exploit this. A biased random basis set of the Zombie states was generated using the biasing regime previously described. Equation (34) was used to ensure that the “dead” and “alive” coefficients were normalized, values θj were randomly generated using a normal distribution centered at 12π for electrons “alive” in the restricted Hartree–Fock Slater determinant [j = 1…6] and centered around 0 for the “dead” electrons. The width, σ, of the normal distributions used to bias each electron is summarized in Table I, and these values are used to create the biased bases in all subsequent results. There are six electrons in the neutral Li2 molecule and four core inactive electrons that are set as always “alive.” Figure 4 schematically shows how this biasing method is applied: the first two spatial orbitals (first four spin orbitals) are set to always be occupied with the final three spatial orbitals (last six spin orbitals) having fractional occupancy.

TABLE I.

The center μ and width σ of the normal distributions used to generate θj for the jth electron of each Zombie state when creating a biased basis.

Figure
j-th electronμ\2πσ\2π
1, 2, 3, 4 0.25 
5, 6 0.25 0.175 
7, 8 0.351 
9, 10 0.120 
Figure
j-th electronμ\2πσ\2π
1, 2, 3, 4 0.25 
5, 6 0.25 0.175 
7, 8 0.351 
9, 10 0.120 
FIG. 4.

Occupational probability plotted against orbital number for five spin orbitals. The first four spin orbitals are set exactly as occupied, and the rest are set using the Gaussian described in Table I. The average probability is shown by a cross, and one standard deviation either side is the bar.

FIG. 4.

Occupational probability plotted against orbital number for five spin orbitals. The first four spin orbitals are set exactly as occupied, and the rest are set using the Gaussian described in Table I. The average probability is shown by a cross, and one standard deviation either side is the bar.

Close modal

Figure 5 also demonstrates that ZS propagation can yield the ground state of the Li2 anion. A 64 biased Zombie state basis was generated, but the first basis function was set as the seven electron open-shell restricted Hartree–Fock Zombie state and the remaining 63 Zombie states generated using the same biasing method as before. Similarly, the first Zombie state can be set to the six-electron restricted Hartree–Fock Zombie state. It can be seen in Fig. 5 that depending on the choice of the initial condition, the same ZS propagation can lead the ground state of neutral molecule and its anion. Figure 6 shows the result of the propagation using the basis of 10, 30, and 50 random basis ZS functions together with their improved energies obtained by cleaning. The first basis function in each basis set is the six electron restricted Hartree–Fock determinant to ensure each evolution starts at the same energy. As expected, the energies for these smaller biased bases are not as accurate as obtained using the “complete” 64 basis function Hamiltonian. However, the biasing method is considerably more accurate when a directly comparing random and biased bases of the same size.

FIG. 5.

Imaginary time propagation for two biased bases each made up of 63 biased Zombie state basis functions and then a Zombie state set either to the six electron restricted Hartree–Fock determinant (orange bottom line) or the seven electron open-shell restricted Hartree–Fock determinant (red top line). The basis with six electron determinant evolves to the Li2 neutral ground state, and the basis with seven electrons evolves to the higher Li2 anion ground state. The corresponding energies, obtained by diagonalizing the complete Slater determinant basis, are also shown as horizontal lines. The details of how each electron is biased are summarized in Table I.

FIG. 5.

Imaginary time propagation for two biased bases each made up of 63 biased Zombie state basis functions and then a Zombie state set either to the six electron restricted Hartree–Fock determinant (orange bottom line) or the seven electron open-shell restricted Hartree–Fock determinant (red top line). The basis with six electron determinant evolves to the Li2 neutral ground state, and the basis with seven electrons evolves to the higher Li2 anion ground state. The corresponding energies, obtained by diagonalizing the complete Slater determinant basis, are also shown as horizontal lines. The details of how each electron is biased are summarized in Table I.

Close modal
FIG. 6.

Imaginary time propagation for basis sets containing 50, 30, and 10 basis functions, the first being the six electron restricted Hartree–Fock determinant and the rest being biased Zombie state basis functions. The details of how each electron is biased are summarized in Table I. The eigenvalue for the Li2 ground state obtained by diagonalizing the Slater determinant basis is given, as the solid horizontal line, for comparison. The cleaned energies, found by E(6e)/N(6e), are shown as dashed lines.

FIG. 6.

Imaginary time propagation for basis sets containing 50, 30, and 10 basis functions, the first being the six electron restricted Hartree–Fock determinant and the rest being biased Zombie state basis functions. The details of how each electron is biased are summarized in Table I. The eigenvalue for the Li2 ground state obtained by diagonalizing the Slater determinant basis is given, as the solid horizontal line, for comparison. The cleaned energies, found by E(6e)/N(6e), are shown as dashed lines.

Close modal

Figure 7 shows the distribution of energy and norm over the Fock spaces with me electrons for the case of 30 BFs. Almost all norm is in the 6e Fock space, but even in this case, cleaning improves the energy and brings it closer to the full CI result. Figure 5 shows that a basis of 64 randomly selected basis functions is complete and yields the exact result even without cleaning.

FIG. 7.

Energy (top) and norm (bottom) distribution for each number of electrons for a basis set of 30 biased Zombie states, the first being set as the six electron restricted Hartree–Fock determinant. The majority of the final wave function belongs to the Fock space on ne = 6 electrons.

FIG. 7.

Energy (top) and norm (bottom) distribution for each number of electrons for a basis set of 30 biased Zombie states, the first being set as the six electron restricted Hartree–Fock determinant. The majority of the final wave function belongs to the Fock space on ne = 6 electrons.

Close modal

Gram–Schmidt orthogonalization was then applied to the complete random basis and is shown in Fig. 8. Four separate states were computed, all of which were initially set to arbitrary vectors. The Gram–Schmidt process was applied and then each state was propagated in imaginary time with Gram–Schmidt orthogonalization applied after every time step. It can be seen in Fig. 8 that the first state is the neutral ground state, the second and third states are the degenerate anion ground states (degenerate as Ms can be ±1/2), and the fourth is the first neutral excited state for Li2. All of these values match the corresponding eigenvalues found by diagonalizing the complete Slater determinant basis.

FIG. 8.

Imaginary time propagation using Gram–Schmidt orthogonalization and the complete random basis for Li2. The final energies of each state correspond accordingly: State 1 is the neutral ground state, states 2 and 3 are the degenerate anion states, and state 4 is the first neutral excited state. These values equal the eigenvalues obtained by diagonalizing the complete Slater determinant basis, which have been shown with the dashed lines.

FIG. 8.

Imaginary time propagation using Gram–Schmidt orthogonalization and the complete random basis for Li2. The final energies of each state correspond accordingly: State 1 is the neutral ground state, states 2 and 3 are the degenerate anion states, and state 4 is the first neutral excited state. These values equal the eigenvalues obtained by diagonalizing the complete Slater determinant basis, which have been shown with the dashed lines.

Close modal

The same Gram–Schmidt process was then applied to four states using a 64 Zombie state biased Hamiltonian. The imaginary time propagation of these four states can be seen in Fig. 9. As with complete random basis propagating in imaginary time with Gram–Schmidt produces the lowest neutral state, the degenerate anion ground states, and the first excited energy, which are also shown on both figures as dashed lines.

FIG. 9.

Imaginary time propagation using Gram –chmidt orthogonalization and the 64 Zombie state biased basis functions for Li2. The final energies of each state correspond accordingly: State 1 is the neutral ground state, states 2 and 3 are the anion ground state, and state 4 is the first excited state. These values equal the eigenvalues obtained by diagonalizing the complete Slater determinant basis, which have been shown with the dashed lines.

FIG. 9.

Imaginary time propagation using Gram –chmidt orthogonalization and the 64 Zombie state biased basis functions for Li2. The final energies of each state correspond accordingly: State 1 is the neutral ground state, states 2 and 3 are the anion ground state, and state 4 is the first excited state. These values equal the eigenvalues obtained by diagonalizing the complete Slater determinant basis, which have been shown with the dashed lines.

Close modal

The findings of this paper are broad but substantially advance Zombie states as a potential method for simulating electronic structure and dynamics. While previously shown to be numerically consistent, the alternative formalism of Zombie states, presented here, builds from the vacuum state using creation and annihilation operators to create Zombie states. This forms a clearer link between Zombie states and more traditional electronic structure theory. This new formulation also leads to a stricter normalization condition, which not only demonstrated how ideas from bosonic systems can be translated to fermionic systems but was key in the implementation of the biasing method used. A key facet of a practicable modern methods is efficient algorithms for the calculation of the Hamiltonian and other operators. In this paper, we have presented algorithms that greatly improve on the naïve algorithms for the number operator, spin operators (Ŝz and Ŝ2), and the two-electron Hamiltonian. A scaled algorithm for the two-electron Hamiltonian is presented and reduces computational cost significantly. Computationally inexpensive calculation of important system properties is obviously very important for any subsequent work. However, this also sets up an easily replicable method for reducing the scaling of any algorithms that may be developed in subsequent work. All of the algorithms follow a similar method using the action of creation and annihilation operators to obtain recursion relations, which lead to lower-scaling algorithms.

We have also shown that using imaginary time propagation on a system of Zombie states, it is possible to find the lowest-lying state of that system. A complete basis of Slater determinants and a basis of complete random Zombie states were shown to give the exact same ground state energies for Li2. Imaginary time propagation removes the need for long real-time evolution and Fourier transforms. Reducing the size of the random basis resulted in the ground state energy obtained no longer matching the benchmark set by the complete Slater determinant basis. However, this problem is overcome by the use of the biasing method, where we use the ideas of CASSCF30 (or CASCI) by splitting orbitals into inactive, active, and virtual types and then biasing individual orbitals in the active space. The lowest four spin orbitals were set as core orbitals and so always “alive,” leaving two electrons to fill the remaining six spin orbitals, such that there were 26 = 64 Slater Determinants in a complete basis set for this active space and any more basis functions would lead to linear dependencies and a singularities in the overlap matrix. However, with 64 biased Zombie state basis functions, it was possible to produce, with imaginary time evolution, energy levels that matched the 16 times larger Slater determinant basis. We also showed that with a single Zombie state changed to a specified restricted Hartree–Fock determinant, it was possible to find a specific energy level with imaginary time propagation. It should be noted that the result of converged imaginary time propagation does not give the correct number of electrons or spin. However, applying cleaning, by construction, does recover this information. Reducing the number of biased Zombie states gives less accurate energy levels; however, when compared to imaginary time propagation for random bases of the same size, the biased bases are considerably better. In some respects, the accuracy of calculations with 64 basis functions is not a great surprise from a chemistry viewpoint, since the lowest four electrons are in 1s orbitals and, therefore, tightly bound to the lithium nuclei. It would be energetically very costly to excite these electrons (compared to the 2s electrons), and they are therefore expected to remain core 1s in any low-lying states. Clearly, the use of Zombie states is not giving new scientific insight into the electronic structure of Li2, but we have demonstrated that ZS can reproduce results with a reduced basis in an efficient manner.

It is worth noting that 64 basis functions is less than the 210 configurations possible with six electrons in ten spin orbitals but far more than the 15 possible configurations possible when the first four electrons are fixed, and the remaining two can be arranged across six spin orbitals. Thus, introduction of Zombie states with non-integer occupations of orbitals for a small system has not decreased the required number of basis functions but, in fact, increased them. This is not surprising as the basis of Zombie States describes multi-electronic states with all possible numbers of electrons from no electrons at all to electrons occupying all spin orbitals. For systems with small active space, such as Li2, the price of not sticking to the right number of electrons outweighs the benefit of bringing in other excited electronic configurations, which may be disregarded by active space methods. However, as can be seen from the Fig. 6, a small ZS basis accounts for some correlation energy. For larger systems, where the active space is large, this may be an affordable way of taking electron correlation into account. The exponential scaling of required configurations is a common problem in electronic structure calculations, but ZS should reduce this, with adequate biasing, with each ZS containing many configurations. As ZSs are made up of superpositions of many electronic configurations, we expect that a manageable number of them will be required to describe larger systems. Therefore, future work will focus on applying Zombie states to much larger systems that traditionally require extremely large basis sets, such as N2,6 to obtain FCI values while using significantly smaller ZS basis sets. We envisage that this work will require the application of sampling techniques previously developed for bosonic systems,33 which when applied to Zombie states should aid in converging results and address the possibility of choosing a linear combination of Zombie states such that the overall wavefunction is an eigenstate of the number operator.

See the supplementary material for details of the code racing for all of the algorithms presented. There is also detail of the code used to generate the electron integrals using PySCF.

The authors would like to thank George Booth for discussions on PySCF and Dmitry Makhov for help with Molpro. The support of the EPSRC via Grant No. EP/P021123/1 is acknowledged. T.J.H.H. also acknowledges a Royal Society University Research Fellowship (Grant No. URF\R1\201502).

The authors have no conflicts to disclose.

O.A.B. and T.J.H.H. contributed equally to this work.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Although the observable properties of a Zombie state can be formally computed exactly by iterative application of creation and annihilation operators followed by computation of the overlap of the resultant states,16 this usually is not an efficient calculation. Here, we present algorithms for the commonly used properties of electronic states, such as the number of electrons, and their spin properties.

If we consider the action of the number operator,

N̂=m=1Mn̂m=m=1Mb̂mb̂m,
(A1)

where b̂m is the annihilation operator and b̂m is the creation operator for Zombie state m, we see

n̂mζ(b)=a11(b)a12(b)a1(m1)(b)a1m(b)a1M(b)a01(b)a02(b)a0(m1)(b)0a0M(b)
(A2)

such that n̂m effectively “deletes” the coefficient a0m(b). Note that the sign switching cancels itself out. It, therefore, appears that acting with the number operator N̂ is O(M) as the summation is over M terms.

A comparison of Eqs. (4) and (A1) suggests that computing the number of electrons represented by any Zombie state, i.e., ζ(a)N̂ζ(b), would be O(M2) as computing the overlap of two Zombie states is O(M).

However, it is possible to construct an algorithm, which is only O(M) by careful consideration of the summation and the product by adapting an algorithm from Ref. 34. Mathematically, calculation of the action of the number operator is

ζ(a)N̂ζ(b)=l=1Mζ(a)n̂lζ(b)=l=1Mj=1l1a1j(a)*a1j(b)+a0j(a)*a0j(b)a1l(a)*a1l(b)×j=l+1Ma1j(a)*a1j(b)+a0j(a)*a0j(b).
(A3)

We can define

di=a1i(a)*a1i(b)
(A4)

and use Eq. (46c) to define

gi=i=1mfi.
(A5)

This can be used along with Eq. (47b) to find the following recursion relations:

g1=f1,
(A6a)
gl=gl1fl,l=2,3,,M,
(A6b)
hM=fM,
(A6c)
hm=hl+1fl,l=M1,M2,,1.
(A6d)

Computation of Eqs. (A4) and (46c) is clearly O(M) and using the recursion relations and so is computation of {gl} and {hm}. Inserting the aforementioned equations into Eq. (A3),

ζ(a)N̂ζ(b)=l=1Mgl1dlhl+1,
(A7)

which is a summation that can also be computed in O(M) steps. When applied, this new algorithm is significantly faster than the previous method; increasing the number of orbitals by a factor of 10, the new algorithm becomes roughly another 10 times faster than the old one as to be expected from the scaling arguments above. Full details of this code racing can be found in Table 5 of the supplementary material.

Since Zombie states are not usually eigenstates of the number operator, unlike most Slater determinants, the uncertainty in the in number of electrons can be defined as a standard deviation,

σN=N̂2N̂.
(A8)

For N̂2,

ζ(a)N̂2ζ(b)=k=1Mζ(a)N̂n̂kζ(b).
(A9)

Computing n̂kζ(b) is trivial (simply delete the “dead” coefficient for orbital k). We can therefore define ζk(b)=n̂kζ(b) and compute

ζ(a)N̂2ζ(b)=k=1Mζ(a)N̂ζk(b)
(A10)

using the known algorithm for computing the number operator, so this is O(M2) instead of O(M3).

The general idea presented trivially extends to the “Ghost” operator,

Ĝ=m=1Mŝmm=1Mb̂mb̂m,
(A11)

which counts how many of the Zombie states are “dead” (unoccupied), as opposed to the number operator that counts how many of them are alive.

1. Ŝz operator

Orbitals are defined with the assumption that any used in the Zombie state calculation are from a restricted calculation such that for each orbital with spin α (spin up, ms = +1/2), there exists an orbital with the same spatial wavefunction but spin component β (spin down, ms = −1/2). For a system with M spin orbitals, there will consequently be K spatial orbitals where K=M/2,KN. We further define that all spin orbitals with up spin have an odd index m = 1, 3, …, M − 1, and all spin orbitals with down spin have an even index m = 2, 4, …, M such that the kth spatial orbital has α spin orbital m = 2k − 1 and β spin orbital m = 2k, k = 1, 2, …, K.

Using this numbering convention, in second quantization notation,

Ŝz=12k=1Kn̂2k1n̂2k
(B1a)
=12m=1M(1)m1n̂m.
(B1b)

The optimized number operator is then adapted introducing a simple sign change rule to calculate Sz in O(M) steps. A comparison of the naïve algorithm, using creation and annihilation operators, and the optimized algorithm, co-opting the optimized number operator equations, can be found in Table 5 of the supplementary material.

2. Faster Ŝz2 computation

Prima facie, computation of Ŝz2 would be O(M3) using Eq. (B1),

Ŝz2=14kM/2jM/2(n̂2k1n̂2k)(n̂2j1n̂2j)=14kMjMn̂kn̂j(1)j+k,
(B2)

so there is O(M) for the summation over k, another over j, and another to evaluate overlap.

However, extending the adapted number operator algorithm used for Ŝz in O(M) calculation to Ŝz2 is straightforward,

ζ(a)Ŝz2ζ(b)=12mMζ(a)Ŝzn̂mζ(b)(1)m,
(B3)

but n̂mζ(b) gives another Zombie state that we call ζ(b,m) such that

ζ(a)Ŝz2ζ(b)=12mMζ(a)Ŝzζ(b,m)(1)m.
(B4)

The bra-ket can be evaluated in O(M) and, therefore, the action of the Ŝz2 operator in O(M2).

3. Total spin

Usually, in electronic structure calculation, one also wants to know the total spin,1 

Ŝ2=Ŝx2+Ŝy2+Ŝz2
(B5a)
=Ŝ+ŜŜz+Ŝz2.
(B5b)

We can apply the faster Ŝz and Ŝz2 algorithms as detailed above. Ŝ+ and Ŝ are raising and lowering operators,

Ŝ+=k=1Nspaŝ+,k,
(B6a)
Ŝ=k=1Nspaŝ,k,
(B6b)
ŝ+,k=b̂2k1b̂2k,
(B6c)
ŝ,k=b̂2kb̂2k1,
(B6d)

where K is the number of spatial orbitals, viz., K = M/2. We are numbering the orbitals 1, 2, …, K such that the kth spatial orbital has an alpha spin orbital number 2k − 1 and a beta spin orbital number 2k.

We then observe the effect of ŝ+,k and ŝ,k on the wavefunction (where m = 2k),

ŝ+,kζ(b)=b̂2k1b̂2ka11(b)a12(b)a1(m1)(b)a1m(b)a1(m+1)(b)a1M(b)a01(b)a02(b)a0(m1)(b)a0m(b)a0(m+1)(b)a0M(b)
(B7a)
=b̂2k1a11(b)a12(b)a1(m1)(b)0a1(m+1)(b)a1M(b)a01(b)a02(b)a0(m1)(b)a1m(b)a0(m+1)(b)a0M(b)
(B7b)
=a11(b)a12(b)a0(m1)(b)0a1(m+1)(b)a1M(b)a01(b)a02(b)0a1m(b)a0(m+1)(b)a0M(b)
(B7c)

and

ŝ,kζ(b)=b̂2kb̂2k1a11(b)a12(b)a1(m1)(b)a1m(b)a1(m+1)(b)a1M(b)a01(b)a02(b)a0(m1)(b)a0m(b)a0(m+1)(b)a0M(b)
(B8a)
=b̂2ka11(b)a12(b)0a1m(b)a1(m+1)(b)a1M(b)a01(b)a02(b)a1(m1)(b)a0m(b)a0(m+1)(b)a0M(b)
(B8b)
=a11(b)a12(b)0a0m(b)a1(m+1)(b)a1M(b)a01(b)a02(b)a1(m1)(b)0a0(m+1)(b)a0M(b).
(B8c)

In both cases, the sign changes cancel each other out exactly. This means that no sign change is required for evaluating Ŝ+ or Ŝ, saving computational cost. However, there are clearly reductions to be made in the scaling that should reduce the time needed to compute ζ(a)ŝ+,lŝ,kζ(b). The action of ŝ,kζ(b) defined by Eq. (B8) can be defined as ŝ,kζ(b)=ζk(c) (where K = M/2 and m = 2l), which has O(M) scaling, and so it is possible to compute

ζ(a)Ŝ+Ŝζ(b)=k=1Kζ(a)Ŝ+ζk(c)=k=1Kl=1Kζ(a)ŝ+,mζk(c)=k=1Kl=1Kj=1l2a1j(a)*a1j(c)+a0j(a)*a0j(c)a0(l1)(a)*a1(l1)(c)a1l(a)*a0l(c)j=l+1Ma1j(a)*a1j(c)+a0j(a)*a0j(c).
(B9)

In a similar manner to the scaled Hamiltonian with Eq. (46c) and the number operator with Eqs. (A5) and (A4), we define

fi=(a1i1(a)*a1i1(c)+a0i1(a)*a0i1(c))(a1i(a)*a1i(c)+a0i(a)*a0i(c)),
(B10a)
gl=i=1lfi,
(B10b)
hl=i=lKfi,
(B10c)
di=a0(i1)(a)*a1(i1)(c)a1i(a)*a0i(c),
(B10d)

which gives

ζ(a)ŝ+lζk(c)=gl2dlhl+1.
(B11)

Overall, this reduces the scaling from O(M3) to O(M2). This new algorithm was then used alongside the original speed improvements giving a large time speedup (see Table 8 in the supplementary material for full details). However, it should be possible to scale the ζ(a)ŝ+ŝζ(b), removing the need to calculate ζ(c). We can set m = 2l and n = 2k and then define

ŝ+,lŝ,kζ(b)=b̂2l1b̂2lb̂2kb̂2k1×a11(b)a12(b)a1m(b)a1n(b)a1M(b)a01(b)a02(b)a0m(b)a0n(b)a0M(b),
(B12)

which can have the following outcomes:

ŝ+,lŝ,kζ(b)=a11(b)a12(b)0a0n(b)a0(m1)(b)0a1M(b)a01(b)a02(b)a1(n1)(b)00a1m(b)a0M(b)if 2k<2l,a11(b)a12(b)a0(m1)(b)00a0n(b)a1M(b)a01(b)a02(b)0a1m(b)a1(n1)(b)0a0M(b)if 2k>2l,a11(b)a12(b)a1(m2)(b)a1(m1)0a1(m+1)(b)a1M(b)a01(b)a02(b)a0(m2)(b)0a0ma0(m+1)(b)a0M(b)if 2k=2l,
(B13)

which gives

ζ(a)ŝ+,lŝ,kζ(b)=j=1k2a1j(a)*a1j(b)+a0j(a)*a0j(b)a0(k1)(a)*a1(k1)(b)a1k(a)*a0k(b)j=k+1l2a1j(a)*a1j(b)+a0j(a)*a0j(b)a1(l1)(a)*a0(l1)(b)a0l(a)*a1l(b)j=l+1Ma1j(a)*a1j(b)+a0j(a)*a0j(b)if 2k<2l,
(B14a)
ζ(a)ŝ+,lŝ,kζ(b)=j=1l2a1j(a)*a1j(b)+a0j(a)*a0j(b)a1(l1)(a)*a0(l1)(b)a0l(a)*a1l(b)j=l+1k2a1j(a)*a1j(b)+a0j(a)*a0j(b)a0(k1)(a)*a1(k1)(b)a1k(a)*a0k(b)j=k+1Ma1j(a)*a1j(b)+a0j(a)*a0j(b)if 2k>2l,
(B14b)
ζ(a)ŝ+,lŝ,kζ(b)=j=1l2a1j(a)*a1j(b)+a0j(a)*a0j(b)a1(l1)(a)*a1(l1)(b)a0l(a)*a0l(b)j=l+1Ma1j(a)*a1j(b)+a0j(a)*a0j(b)if 2k=2l.
(B14c)

We then have to define the following recursion relations:

fi=(a1(i1)(a)*a1(i1)(b)+a0(i1)(a)*a0(i1)(b))(a1(i)(a)*a1(i)(b)+a0(i)(a)*a0(i)(b)),
(B15a)
ci=a0(i1)(a)*a1(i1)(b)a1i(a)*a0i(b),
(B15b)
di=a1(i1)(a)*a0(i1)(b)a0i(a)*a1i(b),
(B15c)
si=a1(i1)(a)*a1(i1)(b)a0i(a)*a0i(b).
(B15d)

These can all trivially be calculated in O(M) steps. We, then, also define

gl=i=1lfi,
(B16a)
hl=i=lKfi,
(B16b)
t(l,p)=i=lpfi.
(B16c)

Equations (B16a) and (B16b) can be calculated recursively in O(M) steps and Eq. (B16c) in O(M2). Therefore,

ζ(a)Ŝ+Ŝζ(b)=k=1K1l=k+1Kgk1ckt(k+1,l1)dlhl+1+l=1K1k=l+1Kgl1dlt(l+1,k1)ckhk+1+l=1Kgl1slhl+1.
(B17)

The naïve total spin algorithm was made up of three separate algorithms Ŝ+Ŝ, Ŝz, and Ŝz2 that scaled O(M3), O(M2), and O(M3), respectively. The scaled algorithm reduced the scaling of each algorithm to O(M2), O(M) and O(M2) respectively. This scaling improvement resulted in Ŝ2 for 1000 orbitals being calculated over 460 times faster than the naïve algorithm; full code racing details are given in Tables 8 and 9 of the supplementary material.

1.
A.
Szabo
and
N. S.
Ostlund
,
Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory
(
Dover Publications
,
New York
,
1989
).
2.
W. M. C.
Foulkes
,
L.
Mitas
,
R. J.
Needs
, and
G.
Rajagopal
,
Rev. Mod. Phys.
73
,
33
(
2001
).
3.
4.
J. B.
Anderson
,
J. Chem. Phys.
63
,
1499
(
1975
).
5.
G. H.
Booth
,
A. J. W.
Thom
, and
A.
Alavi
,
J. Chem. Phys.
131
,
054106
(
2009
).
6.
D.
Cleland
,
G. H.
Booth
, and
A.
Alavi
,
J. Chem. Phys.
132
,
041103
(
2010
).
7.
G. H.
Booth
and
A.
Alavi
,
J. Chem. Phys.
132
,
174104
(
2010
).
8.
D. M.
Cleland
,
G. H.
Booth
, and
A.
Alavi
,
J. Chem. Phys.
134
,
024112
(
2011
).
9.
G. H.
Booth
,
D.
Cleland
,
A. J. W.
Thom
, and
A.
Alavi
,
J. Chem. Phys.
135
,
084104
(
2011
).
10.
J. C.
Greer
,
J. Chem. Phys.
103
,
1821
(
1995
).
11.
L.
Tong
,
M.
Nolan
,
T.
Cheng
, and
J. C.
Greer
,
Comput. Phys. Commun.
131
,
142
(
2000
).
12.
J. P.
Coe
,
D. J.
Taylor
, and
M. J.
Paterson
,
J. Chem. Phys.
137
,
194111
(
2012
).
13.
J. P.
Coe
and
M. J.
Paterson
,
J. Chem. Phys.
137
,
204108
(
2012
).
14.
J. P.
Coe
,
D. J.
Taylor
, and
M. J.
Paterson
,
J. Comput. Chem.
34
,
1083
(
2013
).
15.
J. P.
Coe
and
M. J.
Paterson
,
J. Chem. Phys.
139
,
154103
(
2013
).
16.
D. V.
Shalashilin
,
J. Chem. Phys.
148
,
194109
(
2018
).
17.
S.
Ray
,
P.
Ostmann
,
L.
Simon
,
F.
Grossmann
, and
W. T.
Strunz
,
J. Phys. A: Math. Theor.
49
,
165303
(
2016
).
18.
J. A.
Green
and
D. V.
Shalashilin
,
Phys. Rev. A
100
,
013607
(
2019
).
19.
Y.
Qiao
and
F.
Grossmann
,
Phys. Rev. A
103
,
042209
(
2021
).
20.
J. R.
Klauder
,
Coherent States: Applications in Physics and Mathematical Physics
(
World Scientific
,
Singapore
,
1985
).
21.
J.-P.
Blaizot
,
Quantum Theory of Finite Systems
(
MIT Press
,
Cambridge, MA
,
1986
).
22.
N. S.
Blunt
,
S. D.
Smart
,
G. H.
Booth
, and
A.
Alavi
,
J. Chem. Phys.
143
,
134117
(
2015
).
23.
A. I.
Streltsov
,
O. E.
Alon
, and
L. S.
Cederbaum
,
Phys. Rev. A
81
,
022124
(
2010
).
24.
O.
Bramley
,
C.
Symonds
, and
D. V.
Shalashilin
,
J. Chem. Phys.
151
,
064103
(
2019
).
25.
D.
Huber
and
E. J.
Heller
,
J. Chem. Phys.
89
,
4752
(
1988
).
26.
D. V.
Shalashilin
and
M. S.
Child
,
Chem. Phys.
304
,
103
(
2004
), part of Special Issue: Towards Multidimensional Quantum Reaction Dynamics.
27.

Providing that the starting wavefunction contains contributions from wavefunctions of the same symmetry and number of electrons as the lowest-energy state.

28.
C. A.
Jiménez-Hoyos
,
T. M.
Henderson
,
T.
Tsuchimochi
, and
G. E.
Scuseria
,
J. Chem. Phys.
136
,
164109
(
2012
).
29.
G. E.
Scuseria
,
C. A.
Jiménez-Hoyos
,
T. M.
Henderson
,
K.
Samanta
, and
J. K.
Ellis
,
J. Chem. Phys.
135
,
124108
(
2011
).
30.
G.
Li Manni
,
S. D.
Smart
, and
A.
Alavi
,
J. Chem. Theory Comput.
12
,
1245
(
2016
).
31.
Q.
Sun
 et al, “
PySCF: The Python-based simulations of chemistry framework
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1340
(
2017
).
32.
O.
Bramley
(
2022
). “
Zombie states code
,” GitHub. https://github.com/OBramley/Zombie.
33.
C.
Symonds
,
J. A.
Kattirtzi
, and
D. V.
Shalashilin
,
J. Chem. Phys.
148
,
184113
(
2018
).
34.
T. J. H.
Hele
, “
An electronically non-adiabatic generalization of ring polymer molecular dynamics
,” MChem thesis,
Exeter College, University of Oxford
,
2011
.

Supplementary Material