A deterministic alternative to the full configuration interaction quantum Monte Carlo method

The scope of traditional approaches to full configuration interaction (FCI) has been limited to simple diatomic molecules,^1,2 and there has been little progress in diagonalizing spaces much larger than a billion determinants.^3–5 However, recent progress in alternative approaches to FCI problems has increased the scope of FCI beyond simple diatomic molecules. Two techniques in particular have been important in this progress, full configuration quantum Monte Carlo (FCIQMC)⁶ and density matrix renormalization group (DMRG).^7–9 Both algorithms provide unique advantages, with DMRG being the definitive method for systems in which one can identify degrees of freedom with low levels of entanglement,^10,11 and FCIQMC showing promise for molecules and extended systems in two or more dimensions.^12,13 The success of DMRG and FCIQMC in quantum chemistry is highlighted by their recent applications to unprecedented large-size determinant spaces while also achieving chemical accuracy.^12,14–22

The FCIQMC method is a useful technique with a few limitations which include biased sampling and comparatively computationally expensive simulations. The biased sampling is a result of the initiator approximation²³ generally employed in FCIQMC calculations, which additionally limits the space in which determinants can be sampled. We will denote the initiator algorithm as iFCIQMC. The initiator approximation can cause systematic errors in iFCIQMC calculations that can be unexpectedly large, as evidenced by comparison to variational DMRG results.²⁴ The need for this approximation is related to the Monte Carlo sampling and not necessarily related to power of the technique which we suggest is finding important determinants.

Here we suggest an alternative to stochastic sampling in favor of a completely deterministic version of the FCIQMC technique that efficiently samples the determinant space. We denote this deterministic algorithm by Adaptive Sampling CI (ASCI). ASCI does not have the same problem with quantum entanglement growth that DMRG has for 2 and 3 dimensional systems, and it also appears faster than FCIQMC in computational efficiency. Our approach also surpasses other techniques for finding energetically important determinant subspaces. The majority of traditional CI methods encode relevant physical degrees of freedom based on excitation levels from a reference determinant.²⁵ These excitation-based methods can also suffer from inaccuracy as they miss important parts of determinant space. Within the CI framework several promising ways to circumvent this problem have been suggested that focus on selected CI approaches where one selects relevant determinants based on different criteria.^26–41 After presenting the ASCI method, we establish a connection between FCIQMC and various selected CI approaches. We then apply the ASCI method to Cr₂, a classic hard problem for many computational electronic structure methods.²⁴ Finally, we demonstrate that excited states are straightforward to calculate with ASCI. We note that the calculation of excited states within stochastic FCIQMC is possible, but requires specialized techniques,⁴² or stochastic orthogonalization between walker sets, which is quite different from the method described here.

In the initial development of FCIQMC, one of the original improvements on the method was to take part of the projection step and make it deterministic.²⁰ In this work, we go further and develop a completely deterministic algorithm. Our approach here is to find important determinants in the same manner as FCIQMC, i.e., to sample determinants based on the absolute value of the ground state wave function amplitudes. The FCIQMC technique was originally presented as a projector method in imaginary time and we use this approach to motivate our method. We start by expanding a wave function in the space of determinants,

and the propagator in imaginary time,

which has an asymptotic solution of a stationary state with $dCidτ=0$⁠. In these equations τ is the imaginary time, i is an index for determinants, D_i are determinants, C_i is the amplitude associated with determinant D_i, and H_ij, is the Hamiltonian matrix element 〈D_i|H|D_j〉. In FCIQMC, the parameter E is a free parameter that controls the population of walkers. Here we will consider E to be the ground state energy or our best approximation thereof. The power of FCIQMC is that it ignores the unimportant parts of determinant space, finds important determinants, and samples them according to their amplitudes. For a stationary state we can solve for the individual coefficients as $Ci=−∑i≠jHijCj(Hii−E)$⁠.

The RHS of this equation captures all aspects of the FCIQMC algorithm. The H_ij in the numerator corresponds to the spawning step (transition moves between determinants), the sign of H_ij in the numerator and specifically summing over positive and negative terms corresponds to the annihilation step (cancelation of positive and negative walkers), and the denominator corresponds to the death/cloning step (adding or removing walkers from the simulation). The key to turning FCIQMC into a non-stochastic algorithm is to remove the stochastic sampling and replace it with a deterministic ranking of important determinants. In both approaches the Hamiltonian matrix elements and a wave function are needed. In FCIQMC the wave function is represented by the distribution of walkers at any given step, whereas for our deterministic algorithm, we use an approximate wave function at each iteration as follows:

In this equation, we have labeled the coefficients of the initial wave function as C⁰, and the output coefficients as C¹. The C⁰ and C¹ coefficients correspond to different determinant spaces, where D¹ is the determinant space that includes everything in D⁰ plus all single and double excitations from the determinants in D⁰. Thus for any good approximation to the ground state wave function we can use Eq. (3) to determine the importance of determinants in a much larger space that what is initially included in D⁰. This C¹ coefficients are essentially a first-order perturbation estimate for CI coefficients in the Epstein-Nesbet perturbation theory.³¹ 

The ASCI algorithm is defined by two determinant subspaces: a core space of size cdets and a target space of size tdets. The core space determines the number of terms j we include in the sum in Eq. (3). This is to say we select cdet determinants with the largest $Cj0$ and consider only those to be non-zero in the sum of Eq. (3). The space to be searched is the set of single and double excitations of the core set of determinants. Since our objective is to find the determinants with the largest amplitudes, we in general only need to search determinants connected to those with large amplitudes. We illustrate the use of this approximation with numerical tests in the Discussion section. The target space contains the top tdet determinants, as determined from the ranking, and is the rank of the matrix diagonalized in each iteration.

Initialize size of the core space, D^core, to size cdets and size of the target space, D^target, to size tdets. Set the starting wave function to the Hartree Fock wave function (C₀ = 1, C_i>0 = 0; E = E_HF).

1. Evaluate the perturbed wave function amplitudes over all the single and double substitutions, D^SD, from the core space, D^core,

   $Ai=∑j≠icoreHijCjHii−E.$

   (4)

2. From D^core and D^SD, use the corresponding coefficients {C_i} and {A_i} to select the tdet largest absolute values to define the new target space, D^target.

3. Form and diagonalize H in the target space.

4. The lowest eigenvalue is the new energy E. The largest cdet amplitudes, by magnitude, define the new core space. If the energy is not converged, return to step (1). For simulations presented here, we halted the simulations when the energy converged to 10⁻⁵ Ha. We have included a convergence plot of Cr₂ for (tdets = 1 000 000, cdets = 5000) in Fig. 1.

The determinants found at the end of the simulation will in general be the most important determinants for the ground state wave function. Unlike FCIQMC this technique has no population control bias, initiator bias, and no sign problem. The technique described here provides an inherently variational energy. However, it is possible to extend the accuracy of the technique with perturbation theory, which comes at the expense of the energy no longer being variational.⁴³ We use the following second order correction to the energy:

We note that other perturbative corrections can be considered.⁴³ For the largest systems considered here, only a modest amount of memory is needed, and all parts except the diagonalization step can be trivially parallelized.

When performing an ASCI calculation the first few steps involve the exploration of higher order excitations. If the starting wave function is a single determinant, then each step increases the maximum number of excitations that have been explored by 2. The coefficients A, as calculated in step (1), can span all single and double excitations from the current wave function and is generally very large. Truncating the coefficients that are calculated for A in the main part of the self-consistent loop allows us to consider applying our algorithm to large systems. To maintain size consistency, the value of tdets and cdets will have to grow with system size. For the systems considered in this work, it was possible to converge the value of cdets, as the energy with respect to this parameter can be extrapolated by running with a few different values.

Apart from the diagonalization step, the most computationally expensive task is determining A and finding its largest values. This can be split among many processors, with the only non-trivial communication occurring during the final aggregation of the final values.

Further improvements are possible by generating the natural orbitals after an initial run and using them to recalculate the electron integrals. The natural orbitals can be generated at any point in the simulation from the current best wave function. The natural orbitals are generally thought to produce compact representations of a wave function.⁴⁴ For a difficult system like Cr₂, we find the use of natural orbitals to be crucial to obtain accurate energies.

Before we present our results we consider numerical tests for our approximation of determining the magnitude of A. In general, the coefficients of the ground state wave function will span many orders of magnitude, and searching over the determinants with the smallest amplitude coefficients is inefficient. Using cdets as a parameter to limit the search to only the important determinants appears to be a well controlled approximation that can be converged for all the system we tested. This approach has similarities to the FCIQMC initiator approximation.²³ 

We demonstrate the accuracy of this approximation by considering the CN molecule in an STO-3G basis at 1.194 Å (6240 determinants). We show that we can find the most important determinants in Hilbert space without having to perform a diagonalization over a larger determinant space. Table I provides a comparison of the most important determinants found with different values of cdets and tdets to those obtained from a full diagonalization over the entire space. For a simulation of 100 core determinants and a target space of 200 determinants (100/200) we found 182 of the top 200 determinants (91%). The remaining 9% of missed determinants were found to be close in amplitude to the determinants that replaced them in the target space. Similar result can be seen for all the simulations presented. These results suggest that for some simulations the search algorithm is not highly dependent on the size of the core space and the small percentage of determinants that are missed by the algorithm are replaced by determinants that are similar in importance. Thus we argue that extremely high accuracy is not needed in determining the ranking order. The errors of this approximation is even further reduced when perturbation corrections are used, which can correct for any important determinants that were missed.

Our main goal in presenting this analysis is to demonstrate that a deterministic method is capable of exploring determinant space in a similar manner to FCIQMC. We did not focus on algorithmic speed and we only provide timings here to suggest an upper bound of what can be expected. For the results presented here we built our own implementation of this algorithm and incorporated tools from various codes^45–47 and electron integrals from various packages.^48,49

For C₂ simulations, the convergence of energies to chemical accuracy was easily achieved, and neither the use of natural orbitals nor perturbation corrections are needed. In comparison with the exact results,⁴ we were able to achieve an accuracy of 1 mHa using a diagonalization no larger than 200 000 determinants. The total computer time for a simulation of this size was less than 2 cpu hours. Results are presented in Table II, using different values of tdets and cdets. We also calculated the binding curve for C₂ with our method and compare to FCI results⁵⁰ in Fig. 2. The FCI results were reported to 0.1 mHa accuracy in this previous benchmark. We compare to a subset of these points and our ASCI+PT2 results agree exactly to this level of accuracy. We also ran a cc-pVTZ simulation at the equilibrium distance, and found the ASCI+PT2 energy to be −75.809 19 Ha in comparison to the DMRG(M = 6000) result of −75.809 285 Ha, a difference of 0.1 mHa.

Fig. 3 shows the convergence of our results for Cr₂. In order to make a comparison with previous studies,²⁴ Cr₂ calculations were carried out with the SV basis set⁵¹ at 1.5 Å with 24 active electrons in 30 orbitals and a frozen core. A compact representation of the wave function in this system is dependent on having a good set of orbitals. As part of our algorithm, we run a preliminary calculation with a target space of 100 000 determinants, after which we calculated the natural orbitals. The resulting natural orbitals were used to recalculate the integrals for the production run. The total energy for our most accurate simulation converged to within 16 mHa of the predicted full CI basis set energy.²⁴ 

A perturbation theory analysis⁴³ was performed using Equation (5), which results in a final energy within 1 mHa of the predicted exact result. The timing for the largest simulation (tdets = 10⁶), including the initial run for calculating the natural orbitals but not including the perturbation correction, requires approximately 7 cpu hours when run on a single core of a 2.40 GHz Intel Xeon processor. We expect that further improvements will speed up the ASCI simulations significantly.

In Fig. 4, we plot the binding curve of Cr₂ for DMRG(M = 1600) and ASCI+PT2 for (24e,30o). The DMRG(M = 1600) calculations are taken from a previous benchmark study.⁵³ We also plot the ASCI+PT2 for the non-frozen core calculation of (48e,42o).⁵² The ASCI results were calculated with (tdets = 300 000, cdets = 30 000). At equilibrium both ASCI results match the extrapolated DMRG results to less than 1 mHa accuracy, as given in a previous reference.²⁴ For the (24e,30o) calculation, ASCI is more accurate than the DMRG(M = 1600). As the bond is stretched, it is clear that ASCI is losing some accuracy with regards to DMRG(M = 1600).

To demonstrate the distribution of the determinants located by ASCI we plot a histogram of excitations from the dominant determinants for Cr₂ in Fig. 5. The ratio between different excitations sectors does not change much in increasing tdet values from 10⁵ to 10⁶. For this range, the quadruple excitations generally make up half the wave function, while higher excitations above the quadruples make up roughly 30% of the wave function.

We are also able to calculate excited states with our technique, as they are obtained automatically within the diagonalization procedure. For a 6-31G* basis of C₂, we compare against previous FCI simulations for the first two excited states.⁵⁴ For a simulation at a distance of 1.25 Å with cdets =10⁴ and tdets = 800 000, we have the following energies for the ground state and the first two excited states (−75.7256, −75.6345, −75.6271) in units of Ha. The exact results are (−75.725 995, −75.636 861, −75.628 883).⁵⁴ Thus, although the accuracy of the excited states is not as good as the ground state, it is still straightforward to achieve chemical accuracy (which is generally defined as 0.0016 Ha or 1 kcal/mol).

We can improve the accuracy of excited states by noting that any eigenstate can be used in Eq. (3). Thus we can find and rank determinants by their importance to individual excited states. The excited state optimization can be done simultaneously with the ground state method, or in a state-by-state bootstrap method. For the simultaneous optimization algorithm, we determine a set of important determinants to retain for each excited state. There is a separate search step for each excited state, but one diagonalization step that combines all retained determinants. In contrast, the bootstrap method would converge each excited state one by one, where at each step determinants would be added in specifically for the targeted excited state. The use of natural orbitals averaged over various excited states may also be used to improve the algorithm.^37,55 A detailed study of the targeted excited state technique will be presented elsewhere.

As mentioned earlier, the ASCI method may also be considered a variant of selected CI techniques^31,32,35,39 in which there is considerable current interest.^{26,41,43,56,57} Our method employs first-order perturbation coefficients for selecting determinants and is thus related to the CI method, CIPSI (Configuration Interaction by Perturbation with multiconfigurational zeroth-order wave functions Selected by Iterative process).³¹ Another related technique is the Λ+SD-CI²⁶ method which uses a one-step energy criteria and a one-step approach together with our Eq. (3) in order to find important determinants. Despite these similarities, none of these algorithms have been pushed to achieve chemical accuracy for strongly correlated systems and do not appear to have been benchmarked against FCIQMC or DMRG. The largest Λ+SD-CI calculations included roughly 50 000 determinants and attained 1–3 mHa accuracy for C₂ 6-31G* (comparable to our results in Table II). As shown in this work, it is easy to go more than an order of magnitude in accuracy using our iterative scheme, without significantly increasing the computational effort. The largest selected CI techniques we are aware of have been extended up to 4 × 10⁶ determinants,⁵⁸ but applied to systems for which no benchmarks exist. Monte Carlo Configuration Interaction method (MCCI), is another related technique that can be used to solve CI problems with random sampling in determinant space, however, it has been benchmarked far less than other techniques.^59,60

For the future development of ASCI and other selected CI techniques, it is important to consider how such methods are different from standard CI methods. The difference is largely due to the construction of the Hamiltonian. Selected CI techniques need unique data structures in order to construct the Hamiltonian efficiently.^58,61 A previous study demonstrated that much larger scale simulations, than what we presented here, is possible for selected CI techniques.⁵⁸ We are currently considering various data structures used previously^58,62 and new structures, to determine the best way to scale up our simulations.

We also note that subsequent to posting of this work to the arxiv,⁶³ another group contacted us about improving the techniques presented here and their work has since been published.⁶⁴ Their results show interesting behavior with regards to how many determinants are needed to converge a calculation with respect to increasing system size. Understanding this behavior will ultimately be important in determining just how far we can push ASCI to simulate larger systems.

We have shown that the underlying dynamics of FCIQMC can be used to generate a deterministic algorithm that can be efficiently used to calculate both ground and excited states of chemical systems. We have applied this technique to a known difficult problem in electronic structure theory, the Cr₂ molecule, and shown that chemical accuracy can be achieved with the cpu power available on any modern computer. Our results suggest that the ASCI method (and selected CI methods in general) should be considered as a state of the art CI method in both accuracy and efficiency. The importance of ASCI and FCIQMC is evident as DMRG and post-DMRG methods are not yet well suited for simulations in two and three dimensions. The ASCI method also distinguishes itself from FCIQMC in that excited states and other properties, such as the 2-RDM, are inherently easy to calculate.^65,66

N.M.T. would like to thank Hitesh Changlani, Cyrus Umrigar, Adam Holmes, Bryan O’Gorman, Bryan Clark, and Jonathan Moussa for useful discussions. This work was supported through the Scientific Discovery through Advanced Computing (SciDAC) program funded by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Basic Energy Sciences. We used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation Grant No. OCI-1053575 and resources of the Oak Ridge Leadership Computing Facility (OLCF) at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC0500OR22725. This work was posted to the Cornell preprint archive on March 8, 2016, with identifier arXiv:1603.02686.