NECI: N-Electron Configuration Interaction with an emphasis on state-of-the-art stochastic methods

NECI started off in the late 1990s as an exact diagonalization code for model quantum dots^1,2 and has evolved into a code to perform stochastic diagonalization of large fermionic systems in finite but large quantum chemical basis sets using the Full Configuration Interaction Quantum Monte Carlo (FCIQMC) algorithm.³ This algorithm samples Slater determinant (i.e., antisymmetrized) Hilbert spaces using signed walkers by propagation of the walkers through stochastic application of the second-quantized Hamiltonian onto the walker population. In philosophy, it is similar to the continuum real-space Diffusion Monte Carlo (DMC) algorithm. However, unlike DMC, no fixed node approximation needs to be applied. Instead, the nodal structure of the wave function, as encoded by the signed coefficients of the sampled Slater determinants (SDs), emerges from the dynamics of the simulation itself. However, being based on an FCI parameterization of the wave function, the FCIQMC method exhibits a steep scaling with the number of electrons and is thus only suited for relatively small chemical systems compared to those accessible to DMC. While the common energy measures in FCIQMC methods, namely, the projected, trial energies (cf. Sec. IV) and the energy “shift,” are not variational, a variational energy can be computed from two parallel FCIQMC calculations either directly (cf. Sec. VI) or via the reduced density matrix (RDM) based energy estimator (cf. Sec. VII).

There are also similarities between the FCIQMC approach and the Auxiliary-Field Quantum Monte Carlo (AFQMC) method,^4–6 both being stochastic projector techniques formulated in second quantized spaces. The latter however works in an over-complete space of non-orthogonal Slater determinants and relies on the phase-less approximation⁷ to eliminate the phase problem associated with the Hubbard–Stratonovich transformation of the Coulomb interaction kernel, the quality of this approximation being reliant on the trial wave function used to constrain the path. The objective of AFQMC is the measurement of observables such as the energy by sampling over the Hubbard–Stratonovich fields. On the other hand, FCIQMC works in a fixed Slater determinant space and relies on walker annihilation to overcome the fermion sign problem. The phase-less approximation renders the AFQMC method polynomial scaling with an uncontrolled approximation, while i-FCIQMC, which is, in principle, an exact method, remains exponential scaling. Finally, FCIQMC provides a direct measure of the CI amplitudes of the many-body wave function expressed in the given orbital basis from which observables can be computed including elements of the reduced density matrices (which do not commute with the Hamiltonian) via pure estimators. Exact symmetry constraints, including total spin, can be incorporated into the formalism.⁸ In this sense, the FCIQMC method is closer in spirit to multi-reference CI methods used in quantum chemistry to study multi-reference problems rather than the AFQMC method.

In its original formulation, the algorithm is guaranteed to converge onto the ground-state wave function in the long imaginary-time propagation limit, provided that a sufficient number of walkers are used. This number is generally found to scale with the Hilbert space size and is a manifestation of the sign-problem in this method, essentially implying an exponential memory cost in order to guarantee stable convergence onto the exact solution. In the subsequent development of the initiator method (i-FCIQMC),⁹ this condition was relaxed to allow for stable simulations at relatively low walker populations, much smaller than the full Hilbert space size, albeit at the cost of a systematically improvable bias. While the initiator adaptation removes the strict need for a minimum walker number, it does not eliminate the exponential scaling of the method such that calculations become more and more challenging with increasing system size. To give an idea of the capabilities of the NECI implementation, estimates for the accessible system sizes are given below. The rate of convergence of the initiator error with the walker number has been found to be slow for large systems. This is a manifestation of size-inconsistency error, which generally plagues linear configuration interaction methods. A very recent development of the adaptive shift method¹⁰ mitigates this error substantially, enabling near-FCI quality results to be obtained for systems as large as benzene.

The development of the semi-stochastic method by Petruzielo et al.¹¹ and its further refinements¹² dramatically reduced the stochastic noise and hence improved the efficiency of the method.

The FCIQMC algorithm, as well as its semi-stochastic and initiator versions, is scalable on large parallel machines, thanks to the fact that walker distribution can be distributed over many processors with a relatively small communication overhead. The methods, however, are not embarrassingly parallel owing to the annihilation step of the algorithm (see also Fig. 1). For this reason, parallelization over very large numbers of processors is a highly non-trivial task, but substantial progress has been made, and here, we show that efficient parallelization up to more than 24 000 central processing unit (CPU) cores can be achieved with the current NECI code.

The FCIQMC method has been generalized to excited states¹³ of the same symmetry as the ground state and to the calculation of pure one- and two-particle reduced density matrices via the “replica-trick”^14–17 (and more recently, three- and four-particle RDMs¹⁸). The availability of RDMs enabled the development of the stochastic complete active space self-consistent field (CASSCF) method^19,20 for treating extremely large active spaces. More recently, a fully spin-adapted formulation of FCIQMC has been implemented based on the Graphical Unitary Group Approach (GUGA),⁸ which overcomes the previous limitations of spin-adaptation, which severely limited the number of open-shell orbitals that could be handled. Other advanced developments of FCIQMC in the NECI code include real-time propagation and application to spectroscopy,²¹ Krylov-space FCIQMC,¹² and the similarity transformed FCIQMC,^22–25 which allows the direct incorporation of Jastrow and similar factors depending on explicit electron–electron variables into the wave function.

A number of stochastic methods have been developed as an extension or variation of the FCIQMC approach. These include density matrix quantum Monte Carlo (DMQMC), which allows the exact thermal density matrix to be sampled at any given temperature and also allows straightforward estimation of general observables, including those which do not commute with the Hamiltonian.^16,26 Applications of DMQMC include providing accurate data for the warm dense electron gas.²⁷ Although not implemented in NECI, DMQMC is available in the HANDE-QMC code.²⁸ 

The FCIQMC method has led to the development of a number of highly efficient deterministic selected CI methods, including the adaptive sampling CI method of Head-Gordon and co-workers²⁹ who also established the connection with the much older configuration interaction by perturbation with multiconfigurational zeroth-order wavefunction selected by iterative process (CIPSI) method of Huron et al.³⁰ but with a modified search procedure, while the heat-bath CI method of Holmes and co-workers³¹ was developed from the heat-bath excitation generation for FCIQMC³² together with an initiator-like criterion to select the connected determinants with extreme efficiency. Later, a sign-problem-free semi-stochastic evaluation of the Epstein–Nesbet perturbation energy was developed by Sharma et al.³³ to compute the missing dynamical correlation energy at the second-order in a memory and CPU efficient manner. Other highly related developments of FCIQMC originate in the numerical analysis literature, including the fast randomized iteration,³⁴ further developments by Greene and co-workers,³⁵ and co-ordinate descent FCI of Wang and co-workers.³⁶ 

Depending on the utilized features, the number of electrons and accessible basis sizes can vary. The i-FCIQMC implementation including the semi-stochastic version is highly scalable and has been successfully applied to Hilbert space sizes of up to 10¹⁰⁸ with 54 electrons.³⁷ Atomic basis sets up to aug-cc-pCV8Z for first-row atoms (1138 spin orbitals) are treatable.³⁸ The reduced density matrices can routinely be calculated for use in accurate stochastic-MCSCF¹⁹ for active spaces containing up to 40 electrons and 38 spatial orbitals.^39,40 Real-time calculations are computationally more demanding but can still be performed for first-row dimers using cc-pVQZ basis sets.²¹ For the similarity transformed FCIQMC method, the limiting factor is not the convergence of the FCIQMC, but the storage of the three-body interaction terms imposing a limit of ∼100 spatial orbitals on the currently available hardware.²⁴ The optimized implementations for the application to lattice model systems, such as the Hubbard⁴¹ (in a real- and momentum space formulation), t − J, and Heisenberg models for a variety of lattice geometries, are implemented in NECI. The applicability of FCIQMC to the Hubbard model strongly depends on the interaction strength U/t. For the very weakly correlated regime U/t < 1, FCIQMC is employable up to 70 lattice sites⁴² using a momentum-space basis. In the interesting, yet most problematic, intermediate interaction strength regime in two dimensions, the transcorrelated (similarity-transformed) FCIQMC is necessary to obtain reliable energies in systems up to 50 sites (at and near half-filling).²³ 

The FCIQMC algorithm as implemented in NECI is based on a sparse representation of the wave function and a stochastic application of the Hamiltonian. We start with the full wave function

with coefficients C_i in a many-body basis $Di$⁠. NECI supports Slater determinants or CSFs as a many-body basis; for simplicity, for now, the usage of determinants is assumed, but the algorithm is analogous to CSFs (see also Sec. XII B 2). The FCIQMC wave function is not normalized. The ground state of a Hamiltonian Ĥ is now obtained by iterative imaginary time-evolution, with the propagator expanded to the first order using a discrete time step Δτ such that

which converges to the ground state of Ĥ for τ → ∞ for $|Δτ|<2W$⁠, where W is the difference between the largest and smallest eigenvalues of Ĥ.⁴³ Here, S(τ) is a diagonal shift applied to Ĥ, which is iteratively updated to match the ground-state energy.

The full wave function is stored in a compressed manner, where only coefficients above a given threshold value C_min are stored. Coefficients smaller than C_min are stochastically rounded. That is, in every iteration, a wave function given by coefficients C_i(τ) is stored such that

where $Ci(τ)=Ci$⁠. This compression is applied in every step of the algorithm that affects the coefficients. The value |C_i(τ)| is referred to as the walker number of the determinant $Di$⁠, so $Di$ is said to have |C_i(τ)| walkers assigned.

Applying the Hamiltonian to this compressed wave function is done by separating it into a diagonal and an off-diagonal part as

and then performed in the three labeled steps (a)–(c). First, in the spawning step, the off-diagonal part is evaluated by stochastically sampling the sum over j and storing the resulting spawned wave function as a separate entity, as described in the flow chart in Fig. 1. Then, in the death step, the diagonal contribution is evaluated deterministically, following a stochastic rounding of the resulting coefficients. This step is performed in-place, since the coefficients of the previous iterations are not required anymore. Finally, the spawned wave function from the off-diagonal part is added in the annihilation step, summing up all contributions from the spawned wave function to each determinant. NECI implements the initiator method too, which labels a class of determinants as initiators, typically those with an associated walker number above a given threshold, and effectively zeroes all matrix elements between non-initiator determinants and determinants with C_i(τ) = 0. The implementation thereof is also sketched in Fig. 1.

In the context of FCIQMC calculations, the core functionality of NECI consists of a highly parallelizable implementation of the initiator FCIQMC method⁹ for both real and complex Hamiltonians. There is both a generic interface for ab initio systems, specialized implementations for the Hubbard and Heisenberg models, and the uniform electron gas. The interface for passing input information on the system to NECI is discussed in Sec. XIV. To enable continuation of calculations at a later point, NECI can write the instantaneous wave function and current parameters—such as the shift value—to disk, saving the current state of the calculation.

The NECI program⁴⁴ itself is written in Fortran and requires extended Fortran 2003 support, which is the default for current Fortran compilers. Parallelization is achieved using the Message Passing Interface (MPI),⁴⁵ and support for MPI 3.0 or newer is required. NECI further requires the BLAS⁴⁶ and LAPACK⁴⁷ linear algebra libraries, which are available in numerous packages. The usage of the HDF5 library⁴⁸ for parallel I/O is supported, but not required. If used, the linked HDF5 library has to be built with Fortran support and for parallel applications. For installation, cmake is required, as well as the fypp Fortran preprocessor.⁴⁹ For pseudo-random number generation, the double precision SIMD oriented fast Mersenne Twister (dSFMT)^50,51 implementation of the Mersenne Twister method⁵² is used. The stable version of the program can be obtained from github at https://github.com/ghb24/NECI_STABLE, licensed under the GNU General Public License 3.0. Some advanced or experimental features are only contained in the development version (for access to the development version, contact the corresponding authors). All features presented here are eventually to be integrated to the stable version. The detailed instructions on the installation can be found in the documentation that is available together with the code.

In the following, various important features of NECI are explained in detail. An overview of excitation generation, a fundamental part of every FCIQMC calculation, is given in Sec. II. Then, the semi-stochastic approach (Sec. III), the estimation of energy and use of trial wave functions (Sec. IV), the recently proposed adaptive shift method to reduce the initiator error (Sec. V) and perturbative corrections to this error (Sec. VI), the sampling of the reduced density matrices, which is crucial for interfacing the FCIQMC method with other algorithms (Sec. VII), the calculation of excited states (Sec. VIII), static response functions (Sec. IX) and the real-time FCIQMC method (Sec. X), the transcorrelated approach (Sec. XI), and the available symmetries, including the total spin conservation utilizing GUGA (Sec. XII), are discussed. Finally, the scalability of NECI is explored (Sec. XIII) and the interfaces for usage with other code are presented (Sec. XIV), in particular for the stochastic-MCSCF method (Sec. XV).

A key component of the FCIQMC algorithm is the sampling of the Hamiltonian matrix elements in the spawning step, where the Hamiltonian is applied stochastically. This requires an efficient algorithm to randomly generate the connected determinants with a known probability p_gen for any given determinant, referred to as excitation generation. This typically means making a symmetry constrained choice of (up to) two occupied orbitals in a determinant and (up to) two orbitals to replace them with, such that the corresponding Hamiltonian matrix element is non-zero. If the spin-adapted functions are used rather than determinants, the connectivity rules change, but the main principles are the same.

The spawning probability for a spawn from a determinant $Di$ to a determinant $Dj$ is, in practice, given by

This means that the purpose of p_gen(j|i) of selecting $Dj$ from $Di$ in the spawning probability p_s is to allow the flexibility in the selection of determinants $Dj$ from $Di$ so that, irrespective of how we choose $Dj$ from $Di$⁠, the rate at which transitions occur is not biased by the selection algorithm. In other words, if a particular determinant $Dj$ is only selected rarely from $Di$ (i.e., with low generation probability), then the acceptance of the move (i.e., the spawning probability) will be with correspondingly high probability (i.e., proportional to the inverse of the generation probability). Conversely, if a determinant $Dj$ is selected with relatively high generation probability from $Di$⁠, then its acceptance probability will be correspondingly low. In other words, from the point of view of the exactness of the FCIQMC algorithm, the precise manner in which excitations are made is immaterial: as long as the probability p_gen(j|i) > 0 when |H_ij| > 0, the algorithm will ensure that transitions from D_i → D_j occur at a rate proportional to |H_ij|, and hence, the walker dynamics converges onto the exact ground-state solution of the Hamiltonian matrix. However, from the point of view of efficiency, different algorithms to generate excitations are by no means equivalent.

That is, events with a very large $|Hij|pgen(j|i)$ can lead to very large spawns and thus endanger the stability of an i-FCIQMC calculation. For time-step optimization, NECI offers a general histogramming method, which determines the time step from a histogram of $|Hij|pgen(j|i)$⁠,⁸ and an optimized special case thereof, which only takes into account the maximal ratio.⁵³ If required, internal weights of the excitation generators such a bias toward double excitations are then optimized in the same fashion to maximize the time step.

However, as a result, the time step and thus the overall efficiency of the simulation are driven by the worst-cases of the $|Hij|pgen(j|i)$ ratio discovered within the explored Hilbert space. Thus, an optimal excitation generator should create excitations with a probability distribution to the Hamiltonian matrix elements such that

This is the optimal probability distribution since then, the acceptance rate is solely determined by the time step.³² 

NECI supports a variety of algorithms to perform excitation generation, with the most notable being the pre-computed heat-bath (PCHB) sampling (a variant of the heat-bath sampling presented in Ref. 32, as described in Subsection 3 of the  Appendix), the on-the-fly Cauchy–Schwartz method⁵⁴ (described in Subsection 2 of the  Appendix), the pre-computed Power–Pitzer method,⁵⁵ and lattice-model excitation generators both for real-space and momentum-space lattices. Additionally, a three-body excitation generator and a uniform excitation generator are available, which are essential for treating systems with the transcorrelated ansatz when including three-body interactions.

As heat-bath excitation generation can have high memory requirements, it might be impractical for some systems. There, the on-the-fly Cauchy–Schwartz method can maintain very good $|Hij|pgen(j|i)$ ratios without significant memory cost, albeit at $O(N)$ computational cost, N being the number of orbitals, and possibly with dynamic load imbalance. The details of the Cauchy–Schwartz excitation generation are discussed in the  Appendix.

In many chemical systems, the wave function is dominated by a relatively small number of determinants. In a stochastic algorithm, the efficiency can be improved substantially by treating these determinants in a partially deterministic manner.

Petruzielo et al. suggested a semi-stochastic algorithm,¹¹ where the FCIQMC projection operator $P^=∑ijPij|DiDj|$ is applied exactly within a small but important subspace, which we call the deterministic space, $D$⁠. Specifically, we write

where

The $P^D$ operator therefore accounts for all spawnings that are both from and to determinants in $D$⁠. The stochastic projection operator, $P^S$⁠, contains all remaining terms. The matrix elements of $P^D$ are calculated and stored in a fixed array, and applied exactly each iteration by a matrix–vector multiplication. The operator $P^S$ is then applied stochastically by the usual FCIQMC spawning algorithm.

The semi-stochastic adaptation requires storing the Hamiltonian matrix within $D$⁠, which we denote $HD$⁠. In NECI, $HD$ is stored in a sparse format, distributed across all processes. To calculate $HD$⁠, we have implemented the fast generation scheme of Li et al.⁵⁶ This approach has allowed us to use deterministic spaces containing up to ∼10⁷ determinants. However, a more typical size of $D$ is between 10⁴ and 10⁵.

Ideally, a deterministic space of a given size (⁠$ND$⁠) should be chosen to contain the determinants with the largest value of |C_i| in the exact FCI wave function. This optimal choice is not possible in practice, but various approaches exist to make an approximate selection. Petruzielo and co-workers suggested using selected configuration interaction (SCI) to make the selection.¹¹ Within NECI, the most common approach is to choose the $ND$ determinants that have the largest weight in the FCIQMC wave function at a given iteration.¹² Therefore, a typical FCIQMC simulation in NECI will be performed until convergence (at some iteration number N_conv.) using the fully stochastic algorithm, at which point the semi-stochastic approach is turned on, selecting the $ND$ most populated determinants in the instantaneous wave function to form $D$⁠. The appropriate parameters (⁠$ND$ and N_conv.) are specified in the NECI input file. NECI supports performing periodic re-evaluation of the $ND$ most populated determinants, updating the deterministic space $D$ with a given frequency.

Using the semi-stochastic adaptation with a moderate deterministic space (on the order of ∼10⁴) can improve the efficiency of FCIQMC by multiple orders of magnitudes. This is particularly true in weakly correlated systems. The semi-stochastic approach can also be used in NECI when sampling the reduced density matrices (RDMs), as described in Sec. VII. Here, contributions to RDMs are included exactly between all pairs of determinants within $D$⁠. It has been shown that this can substantially reduce the error on RDM-based estimators.¹² Using the semi-stochastic adaptation in NECI disables the load-balancing unless a periodic update of $D$ is performed.

The most common energy estimator used in FCIQMC is the reference-based projected estimator,

where |D_Ref⟩ is an appropriate reference determinant (usually the Hartree–Fock determinant). In case $Ψ$ is an eigenstate, this yields the exact energy, but, in general, it is a non-variational estimator. This is the default estimator for the energy and can be obtained with a minimal overhead.

NECI has the option to use the projected estimators based on more accurate trial wave functions, which can significantly reduce statistical error in energy estimates. For this reason, we define a trial subspace $T$⁠, which is spanned by $NT$ determinants. Similarly to the deterministic space, $T$ should ideally be formed from the determinants with the largest contribution in the FCI wave function or some good approximation to these determinants. Projecting Ĥ into $T$ gives us a $NT×NT$ matrix, which we denote $HT$⁠, whose eigenstates can be used as trial wave functions for more accurate energy estimators.

Let us denote an eigenstate of $HT$ by $|ΨT=∑i∈TCiT|Di$ with eigenvalue $ET$⁠. Then, a trial function-based estimator can be defined as

Here, $C$ is the space of all determinants connected to $T$ by a single application of Ĥ (not including those in $T$ ). C_i denotes walker coefficients in the FCIQMC wave function, and V_j is defined within $C$ as

To calculate the estimator E_Trial, we therefore require several large arrays: first, $HT$⁠, which is stored in a sparse format, in the same manner as the deterministic Hamiltonian in the semi-stochastic scheme; second, $ψT$⁠, which must be calculated by the Lanczos or Davidson algorithm; third, V, which is a vector in the entire $C$ space. The number of coefficients to store in $C$ is larger than in $T$ by a significant amount, typically by several orders of magnitude. Indeed, storing $V$ can become the largest memory requirement. Because of this, using trial wave functions is typically more memory intensive in NECI than using the semi-stochastic approach for a given space size. We therefore suggest using a smaller trial space, $T$⁠, compared to the deterministic space, $D$⁠.

Note that the initiator error on E_Trial is not the same as the initiator error on E_Ref. For example, E_Trial becomes exact as $|ΨT$ approaches the FCI wave function. However, for practical trial wave functions, the two energy estimates typically give similar initiator errors for ground-state energies in our experience. An exception occurs for excited states (see Sec. VIII). In this case, the wave function is usually not well approximated by a single reference determinant, and E_Trial with an appropriate $T$ yields a great improvement, both for the statistical and initiator errors.

The initiator criterion⁹ is important in making FCIQMC a practical method allowing us to achieve convergence at a dramatically lower number of walkers than the full FCIQMC.³ However, this approximation introduces a bias in the energy when an insufficient number of walkers are used. This bias can be attributed to the fact that non-initiators are systematically undersampled due to the lack of feedback from their local Hilbert space. To correct this, we can allow each non-initiator determinant $Di$ to have its own local shift S_i(τ) as an appropriate fraction of the full shift S(τ),

The fraction f_i is computed by monitoring which spawns are accepted due to the initiator criterion and accumulating positive weights over the accepted and rejected ones,

These weights $wij$ are derived from perturbation theory⁵⁷ where the first-order contribution of determinant $Di$ to the amplitude of determinant $Dj$ is used as a weight for spawns from $Di$ to $Dj$⁠,

It is worth noting that, regardless of how the weights are chosen, expression (14) guarantees that initiators get the full shift. In addition, as the number of walkers increases, the local Hilbert space of a non-initiator becomes more and more populated, restoring the full method in the large walker limit.

We call the above approach for unbiasing the initiator approximation, the adaptive-shift method.¹⁰ In Fig. 2, the exemplary results (from Ref. 10) from using the adaptive shift method are displayed, comparing the total energies of the butadiene molecule in ANO-L-pVDZ basis (22 electrons in 82 spatial orbitals), obtained with the normal initiator method and the adaptive shift method using three different values of the initiator parameter n_a: 3, 10, and 20. The adaptive shift results are in good agreement with other benchmark values from Density Matrix Renormalization Group (DMRG), coupled-cluster singles doubles triples perturbative quadruples [CCSDT(Q)], and extrapolated HCIPT2. In contrast, the normal initiator method has a bias of over 10 mH. Also note how by using the adaptive shift, the results become, to a large extent, independent of the initiator parameter n_a.

An alternative approach to removing initiator error in NECI is through a perturbative correction.⁶⁰ In the initiator approximation, spawning events from non-initiators to unoccupied determinants are typically discarded. These discarded events make up a significant fraction of all spawning attempts made, which, in turn, accounts for much of the total simulation time. While it is necessary to discard these spawned walkers to prevent disastrous noise from the sign problem,⁶¹ this step is extremely wasteful.

These discarded walkers actually contain significant information, which can be used to greatly increase the accuracy of the initiator FCIQMC approach. Specifically, these walkers may sample up to double excitations from the currently occupied determinants (a similar argument can be used to justify the above adaptive shift approach). In analogy with a comparable approach taken in selected CI methods, these discarded walkers can be used to sample a second-order correction to the energy from Epstein–Nesbet perturbation theory.

The correction is calculated by

Here, Δτ is the time step, E₀ is the i-FCIQMC estimate of the energy, and $Sir$ is the total spawned weight onto the determinant |D_i⟩ in replica r (the replica approach will be discussed in more detail in Sec. VII). This correction requires that two replica FCIQMC simulations are being performed simultaneously to avoid biases in this estimator. The summation here is performed over all spawning attempts, which are discarded on both replicas simultaneously.

This must only be applied to correct the variational energy estimator from i-FCIQMC. Such variational energies in NECI can be calculated either directly^62,63 or from the two-body reduced density matrices, which may be sampled in FCIQMC.

This perturbative correction is essentially free to accumulate, since all spawned walkers contributing to Eq. (16) are created regardless. The only significant extra cost comes from the requirement to perform two replica simulations. However, for large systems, the noise on this correction can become significant, which necessitates further running time to reduce statistical errors.

This correction has proven extremely successful in practice, particularly for weakly correlated systems, where it is typical to see 80%–90% of remaining initiator error removed.^60,62,63

While the total energy is an important quantity to extract from quantum systems, a more complete characterization of a system requires the ability to extract information about other expectation values. If these expectation values are derived from operators that do not commute with the Hamiltonian of the system, then a “projected” estimate of the expectation value akin to Eq. (9) is not possible, and alternatives within FCIQMC are required in order to compute them. This is the case for many key quantities such as nuclear derivatives (forces on atoms), dipole moments, and higher-order electrical moments, as well as other observables such as pair distribution functions.⁶⁴ They all can be obtained via the corresponding n-body reduced density matrix (n-RDM), where n is the rank of the operator in question, which fully characterizes the correlated distribution and coherence of n electrons relative to each other. This information can also be used to calculate quantum information measures, which are not observables but which characterize the entanglement within a system, such as correlation entropies.¹⁵ 

To characterize the strength of coupling between different states under certain operators, e.g., the oscillator strength of optical excitations, as well as obtaining other dynamical information requires computing transition density matrices (tRDMs) between stochastic samples of different states, which can be sampled within FCIQMC using the excited state feature discussed in Sec. VIII.^17,65 Furthermore, the two states considered may not sample eigenstates of the system, but one of them can be a response state of the system, then the resulting tRDMs characterize the response of a system to a perturbation, corresponding to a higher derivative of the energy such as the polarizability of the system, which will be addressed in Sec. IX.⁶⁶ Finally, RDMs can also be used to characterize the expectation value of an effective Hamiltonian in a subspace of a system.^67,68 This effective Hamiltonian can include effects such as electronic correlations coupling the space to a wider external set of states. The plurality of electronic structure methods of this kind, such as explicitly correlated “F12” corrections for basis set incompleteness,^69–71 multi-configurational self-consistent field,^19,20 internally contracted multireference perturbation theories,¹⁸ embedding methods,^72,73 and the Multi-Configuration Pair-Density Functional Theory (MC-PDFT),⁷⁴ further attests the importance of faithful and efficient sampling of RDMs in electronic structure theory.

All expectation values of interest can be derived from contractions with a general reduced density matrix object, defined as

where n denotes the “rank” of the RDM and the choice of the states A and B defines the type of RDM, as described above. In this section, we focus on the sampling of the 2-RDM. This is generally the most common RDM required, as most expectation values of interest are (up to) two-body operators, including the total energy of the system. Furthermore, within FCIQMC, the fact that the rank of the RDM required is then the same as the rank of the Hamiltonian, which is sampled within the stochastic dynamics, leads to a novel algorithm, which ensures that the overhead to compute the 2-RDM is relatively small and manageable.¹⁵ 

Expanding the expression for the 2-RDM in terms of the exact FCI wave function [Eq. (1)], we find

where i, j index the many-electron Slater determinants and k, l, m, n denote single-particle orbitals. We will focus on the case where we are sampling |Ψ_A⟩ = |Ψ_B⟩ = |Ψ₀⟩, the ground state of the system, since the same basic principles are applied to sampling the tRDMs, where the other walker distribution may represent an excited state or a response state, with more details for these cases considered in Refs. 17 and 66. The expectation values derived from these RDMs describe “pure” expectation values to distinguish them from the projective estimate of expectation values given in Eq. (9).

There are some features of the form of Eq. (18) that should be noted. First, the 2-RDM requires the sampled amplitudes on all determinants in the space connected to each other via (up to) double electron substitutions. This means that this expectation value requires a global sampling of connections in the entire Hilbert space, in contrast to the projected energy estimate, which requires only a consideration of the determinant amplitudes that are connected directly through Ĥ to the reference determinant (or small trial wave function, see Sec. IV). Second, it is seen that the pairs of determinants in Eq. (18) are exactly the same as the pairs of determinants connected, in general, through the Hamiltonian operator used to sample the FCIQMC dynamics in Eq. (4), assuming that the matrix element is not zero due to (accidental) symmetry between the determinants. This allows an algorithm to sample the 2-RDM concurrently with the sampling of the Hamiltonian required for the spawning steps between occupied determinant pairs.

A final point to note is that the n-RDM is a non-linear functional of the FCI amplitudes—specifically being a quadratic form. Within the FCIQMC sampling, the C_i amplitudes are stochastic variables represented as walkers [C_i(τ)], which, at any one iteration, are, in general, very different from the true C_i but, when averaged over long times, have an expected mean amplitude, which is the same as (or a very good approximation to) C_i. However, due to this non-linearity in the form of the 2-RDM, the average of the sampled amplitude product is not equal to the product of the average amplitude, $⟨Ci*(τ)Cj(τ)⟩τ≠⟨Ci*(τ)⟩τ⟨Cj(τ)⟩τ$⁠, as it neglects the (co-)variance between the sampled determinant amplitudes. Initial applications of RDM sampling in FCIQMC neglected these correlations in the sampling of the RDMs, which significantly hampered the results, especially for the diagonal elements of the RDMs.⁶⁹ The result is that even if each determinant were correctly sampled on average, the stochastic error in the sampling would manifest as a systematic error in the RDMs, and thus, only give correct results in the large walker limit, but not the large sampling limit, even if the wave function were correctly resolved.

The resolution to this problem came via the “replica trick,”^15,16 which changes the quadratic RDM functional into a bilinear one.¹⁴ This formally removes the systematic error in the RDM sampling, at the expense of requiring a second walker distribution. The premise is to ensure that these two walker distributions are entirely independent and propagated in parallel, sampling the same (in this instance ground-state) distribution. This ensures an unbiased sampling of the desired RDM, by ensuring that each RDM contribution is derived from the product of an uncorrelated amplitude from each replica walker distribution. The sampling algorithm then proceeds by ensuring that during the spawning step, the current amplitudes are packaged and communicated along with any spawned walkers. During the annihilation stage, these amplitudes are then multiplied by the amplitude on the child determinant from the other replica distribution, and this product then contributes to all n-RDMs, which are accumulated and equal to the rank of the excitation or higher. In this way, the efficient and parallel annihilation algorithm is used to avoid the latency of additional communication operations, with the necessary packaging of the amplitude and specification of the parent determinant along with each spawned walker being the only additional overhead. The NECI implementation allows for up to 20 replicas to be run, which exceeds any needs arising in the context of RDM calculation.

Full details about the ground-state 2-RDM sampling algorithm can be found in Ref. 15; however, we mention a few salient additional details here. The RDMs are stored in fully distributed and sparse data structures, allowing the accumulation of RDMs for very large numbers of orbitals. The sampling of the RDMs is also not inherently Hermitian. While the sampling within FCIQMC obeys detailed balance, the flux of walkers spawned from $Di→Dj$ is only equal to the reverse flux on average, and therefore, the stochastic noise ensures that the swapping of the two states does not give identical accumulated RDM amplitudes for finite sampling (note that for transition RDMs, this is not expected, with more details in Ref. 17). Similarly, the states sampled in FCIQMC are not normalized, and therefore, neither are the sampled RDMs. Both of these aspects are addressed at the end of the calculation, where the RDMs are explicitly made Hermitian via averaging appropriate entries, and the normalization is constrained by ensuring that the trace of the RDMs gives the appropriate number of electrons.¹⁵ 

The dominant cost of RDM sampling in large systems comes from the sampling of elements defined by pairs of creation and annihilation operators with the same orbital index. These correspond to tuples of occupied orbitals common to both |D_i⟩ and |D_j⟩ states. We term these contributions promotions, as they contribute to a rank of a RDM greater than the excitation level between |D_i⟩ and |D_j⟩. For instance, single excitation spawning events need to contribute to all N − 1 elements of the 2-RDM corresponding to common occupied orbitals in the two determinants. The most extreme case comes from the “diagonal” contributions to the RDMs, where i = j, which requires N(N − 1)/2 contributions to the 2-RDM to be included where each index defining the RDM element corresponds to the same occupied orbital in the two determinants. To mitigate this cost, these diagonal elements are stored locally on each MPI process and only infrequently accumulated at the end of an RDM “sampling block,” or when the determinant becomes unoccupied, with the amplitude averaged over the sampling block. This substantially reduces the frequency of the $O(N2)$ operations required to sample these promoted contributions from the diagonal of Eq. (18).

Other efficiency boosting modifications to the algorithm, such as the semi-stochastic adaptation¹² (detailed in Sec. III), are also seamlessly integrated with the RDM accumulation. Within the deterministic core space, the RDM contributions are exactly accumulated along with the exact propagation, with the connections from the deterministic to the stochastic spaces sampled in the standard fashion. This combination of RDM sampling with the semi-stochastic algorithm can greatly reduce the stochastic errors in the RDMs by ensuring that contributions from large weighted determinant amplitudes are explicitly and deterministically included. Furthermore, the reference determinant and its direct excitations are also exactly accumulated. This is partly because these are likely important contributions, but principally, if the reference is a Hartree–Fock determinant, then the coupling to its single excitations via the Hamiltonian will be zero due to Brillouin’s theorem. These single excitations will nevertheless contribute to the RDMs and therefore are included explicitly.

The sampling of RDMs with a rank greater than two is also now possible within the FCIQMC algorithm and NECI code. The importance of these quantities is primarily in their use in internally contracted multireference perturbation theories, although a number of other uses for these quantities also exist.¹⁸ These methods allow for the FCIQMC dynamics to only consider an active orbital subspace, hugely reducing both the full Hilbert space of the stochastic dynamics and the required time step, while the accumulation of up to 4-RDMs (or contracted lower-order intermediates for efficiency) allows for a rigorous coupling of the strong correlation in the low-energy active space to the dynamic correlation in the wider “external” space via post-processing of these higher-body RDMs with integrals of the external space. Sampling of higher-body RDMs cannot use the identical algorithm to the 2-RDMs, since it now requires the product of determinant amplitudes separated by up to 4-electron excitations, which are not explicitly sampled via the standard FCIQMC propagation algorithm. To allow for this sampling, we include an additional spawning step per walker of excitations with a rank between three and n, where n is the rank of the highest RDM accumulated. This additional spawning is controlled with a variable stochastic resolution, ensuring that the frequency of these samples is relatively rare to control the cost of sampling of these excitations [approximately only one higher-body spawn for every 10–20 traditional (up to two-body) spawning attempts]. There is no time step associated with these excitations, and every attempt is “successful,” transferring information about higher-body correlations in the system and contributing to these higher-body excitations, but not modifying the distribution of the sampled wave function. However, the dominant cost of sampling these higher-body RDMs is not the sampling events themselves, but rather the promotion of lower-rank excitations to these higher-body intermediates. Nevertheless, the faithful sampling of these higher-body properties has allowed for the stochastic estimate of fully internally contracted perturbation theories in large active spaces, with a similar number of walkers required to sample the 2-RDM in an active space.¹⁸ 

In many applications, besides ground-state energies, the properties of excited states are of interest. If states in different symmetry sectors are targeted, this can be easily achieved by performing separate calculations in each sector, yielding the ground state with a given symmetry. If, however, several eigenstates with the same symmetry are required, then this approach is not sufficient. The FCIQMC method is not inherently limited to ground-state calculations and can employ a Gram–Schmidt orthogonalization technique to calculate a set of orthogonal eigenstates.^13,17 The obtained states will then be the lowest energy states with a given symmetry.

Calculating eigenstates sequentially and orthogonalizing against all previously calculated states carry the problem of only orthogonalizing against a single snapshot of the wave function, which will lead to a biased estimate of the excited states. Instead, calculating all states in parallel and orthogonalizing after each iteration give much better results.

The required modifications to the algorithm are minimal. To calculate a set of m eigenstates, m FCIQMC calculations are run in parallel, with the additional step of performing the instantaneous orthogonalization between the m states, performed at the end of each iteration. The orthogonalization requires $O(m2)$ operations and uses one global communication per state. To run m parallel calculations, the replica feature presented in Sec. VII is used to efficiently sample a number of states in parallel. After each FCIQMC iteration, for each state, the contributions from all states of lower energies are projected out. The update step for the nth wave function $ψn$ is then modified to

with the orthogonalization operator for the nth state,

With this definition of the orthogonalization operator, the ground-state FCIQMC wave function (n = 0) is left unaffected. The first excited state (n = 1) is then orthogonalized against the ground state (using the updated wave functions at τ + Δτ, after annihilation has been performed). The second excited state is orthogonalized against both the ground and first excited states, and so on.

To enforce the FCIQMC wave function discretization, after performing the orthogonalization, all determinants with a coefficient smaller than the minimal threshold (typically 1) are stochastically rounded (either down to 0 or up to 1, in an unbiased manner). This is required to prevent proliferation of very small walkers, which adversely affects the wave function compression.

Response theory is a well-established formalism to calculate molecular properties using quantum chemical methods.^75–78 It is, in general, formulated for a time-dependent field, which allows us to compute both static and dynamic molecular properties. However, it is currently only implemented for a static field within NECI.⁶⁶ 

Calculation of molecular properties using response theory relies on the evaluation of the response vectors that are the first or higher order wave functions of the system in the presence of an external perturbation $V^$⁠. According to Wigner’s “(2n + 1)” rule, response vectors up to order n are required to obtain response properties up to order 2n + 1.⁷⁷ For calculating second-order properties such as dipole polarizability, the first-order response vector, C⁽¹⁾, needs to be obtained along with the zero-order wave function parameter C⁽⁰⁾. While C⁽⁰⁾ uses the original FCIQMC working equation (4), C⁽¹⁾ is updated according to

The response vector is discretized into signed walkers in the same way it is done for C⁽⁰⁾. The dynamics of the response-state walker is simulated according to Eq. (21) using an additional pair of replica, and it works in parallel with the dynamics of the zero-order state. Additional spawning and death steps are devised for the response-state walker dynamics due to the presence of the perturbation, alongside the original spawning and death steps in the dynamics. The dependence of the response state on the zero-order states comes from these two aforementioned additional steps. A Gram–Schmidt orthogonalization is applied to the response-state walker distribution with respect to the zero-order walker distribution at each iteration using the same functionality as described in Sec. VIII. This ensures the orthogonality of the response vectors with respect to all lower-order wave function parameters.

The norm of the response walkers is fixed by the choice of the normalization of the zero-order walkers, and it can, in principle, grow at a much faster rate than the zero-order norm. Therefore, in Eq. (21), we introduce the parameter α to control the norm of the response walkers and to reduce the computational effort expended in simulating their dynamics. We aim at matching the number of response-state walkers [$Nw(1)$] with the number of zero-order walkers [$Nw(0)$] by updating α periodically as

Once the walker number stabilizes, the value of α is kept fixed while accumulating statistics. As α scales the norm of the response vector, it needs to be taken into account while evaluating response properties.

Response properties are then obtained from the transition reduced density matrices (tRDMs) that are stochastically accumulated following Eq. (18). For example, dipole polarizability is obtained from the one-electron tRDMs between the zero- and first-order wave functions as

with the $γp,qy$ being calculated from the two-electron tRDM as

Due to the use of two replica per state while sampling both zero- and first-order states, a statistically independent and unbiased estimator of tRDMs can be constructed in two alternative ways, which are denoted here as “[1]” and “[2].” The perturbation used in the computation of the tRDMs in Eq. (24) is the dipole operator ŷ. The factor $1α$ appears due to the re-scaling of the response vector following Eq. (21).

For the purpose of obtaining spectroscopic data or targeting highly excited states, the calculation of orthogonal sets of eigenstates quickly becomes unfeasible, as to obtain a certain eigenstate, all eigenstates of lower energy with the same symmetry have to be computed as well. Spectral functions and the resulting excitation energies can however be calculated using real-time evolution of the wave function, yielding time-resolved Green’s functions, which contain information on the full spectrum. In addition to the stochastic imaginary time evolution of a wave function using in the calculation of individual states, NECI supports performing real-time and arbitrary complex-time calculations, evolving the wave function alongside a complex time trajectory.²¹ As Green’s functions are quadratic in the coefficients of the wave function and averaging over multiple iterations is not an option when evolving a wave function with a real-time component, running multiple calculations in parallel akin to excited state calculations discussed in Sec. VIII is mandatory, as is running with complex coefficients. The real-time propagation can be used to obtain energy gaps from spectral densities and thus target excited states. In contrast to the direct calculation of excited states, these have not to be calculated one by one and in order of ascending energy, however. In Fig. 3, a simple example for applying both the excited-state search and the real-time evolution to the Beryllium atom in a cc-pVDZ basis set to obtain the singlet–triplet gap of the lowest P-state is given. An issue with running real-time calculations is the difficulty of population control, as the death step is essentially replaced by a rotation in the complex plane. This issue can be mitigated by a rotation of the trajectory, evolving along a trajectory in the complex plane. NECI supports an automated trajectory selection that updates the angle α of the time trajectory in the complex plane to maintain a constant population. The Green’s function obtained in the complex plane can then be used to obtain real-frequency spectral functions using analytic continuation,^79,80 with the analytic continuation being significantly easier and more information being recoverable the closer to the real axis the trajectory is.²¹ As, in contrast to the projector FCIQMC, errors arising from the expansion of the propagator are a concern when running complex-time calculations, NECI uses a second-order Runge–Kutta integrator here to sufficiently reduce the time-step error.

The computational cost of a full CI method usually scales exponentially with respect to the size of the basis set. On the other hand, the low regularity of wave functions (characterized by the electronic cusp⁸³) causes a very slow convergence toward the basis set limit. For calculations aiming at highly accurate results, it is very helpful to speed up such slow convergence.

A Jastrow⁸⁴ ansatz offers a way to factor out the cusp from the wave function,

where $T^=∑i<ju(ri,rj)$ is a symmetric function [u(r_i, r_j) = u(r_j, r_i)] over electron-pairs and $Φ$ is an anti-symmetric many-body function. By including the cusp term |r_i − r_j|/2 in u(r_i, r_j), the regularity of $Φ$ is improved at least by one order over $Ψ$⁠.⁸⁵ We can also include other terms in u(r_i, r_j) to capture as much dynamic correlations as possible. By using variational quantum Monte Carlo (VMC) methods, the pair correlation function u(r_i, r_j) can be obtained for a single determinant $Φ$ (e.g., $ΦHF$⁠) or a linear combination of a small number of determinants (e.g., a small CAS wave function).

The transcorrelated method of Boys and Handy⁸⁶ provides a simple and efficient way to treat the Jastrow ansatz, where the original Schrödinger equation is transformed into a non-Hermitian eigenvalue problem,

The advantage of this form of $T^$ is that the similarity transformation leads to an expansion, which terminates at the second order,

The similarity transformation introduces a novel two-body operator $K^$ and a three-body potential $L^$⁠,

The whole transcorrelated Hamiltonian can be written in the second quantized form as

where $hqp=ϕp|h|ϕq$ and $Vrspq=ϕpϕq|r12−1|ϕrϕs$ are the one- and two-body terms of the molecular Hamiltonian, while $Krspq=ϕpϕq|K^|ϕrϕs$ and $Lstupqr=ϕpϕqϕr|L^|ϕsϕtϕu$ originate from the $K^$ and $L^$ operators.

This transcorrelated method has been investigated by FCIQMC using NECI, as it can essentially speed up the convergence with respect to basis sets. On the other hand, the effective Hamiltonian is non-Hermitian and contains up to three-body potentials. Luo and Alavi explored a transcorrelated approach where only up to two-body potentials are involved.²² The performance on uniform electron gases indicates that this approach could be developed into an efficient FCIQMC method for plane wave basis sets in the future. For general molecular systems, the full transcorrelated Hamiltonian (32) is implemented in NECI, where $T^$ is fixed and treated as an input function, while $Φ$ is sampled by the FCIQMC algorithm. The lack of a lower bound of the energy due to the non-Hermiticity of the similarity transformed Hamiltonian poses a severe problem for variational approaches. However, as a projective technique, FCIQMC does not have an inherent problem sampling the ground-state right eigenvector by repetitive application of the projector (2) and obtaining the corresponding ground-state eigenvalue.

The matrix elements $Krspq$ and $Lstupqr$ are pre-calculated and have to be supplied as an input. The matrix elements of K can be passed combined with the Coulomb integrals, while the matrix elements of L are passed in a separate input file. This treatment is efficient for small atomic and molecular systems, but for large systems, the storage of the L matrix becomes a bottleneck. Here, efficient low rank tensor product expansion of L might in the future make it practical to treat even larger systems. NECI supports storage of L in a dense and a sparse format as well as on-the-fly calculation of $Lstupqr$ from a tensor decomposition. Additionally, major technical changes to the FCIQMC implementation are required for sampling up to triple excitations, which generally leads to reduced time steps. The development of efficient excitation generation, which can alleviate the time-step bottleneck, is the subject of current work.

This method has been tested on the first row atoms,²⁴ which shall serve as an example here. Two different correlation factors obtained by Schmidt and Moskowitz⁸⁷ based on variance-minimization VMC, which contain 7 and 17 terms of polynomial type basis functions, have been employed there. The 7-term factor (SM7) contains mainly electron–electron correlation terms together with some electron–nuclear terms, while the 17-term factor (SM17) uses more terms to describe also the electron–electron–nuclear correlation. For the full CI expansion of $Φ$⁠, three different basis sets, cc-pVDZ, cc-pVTZ, and cc-pVQZ, respectively, have been used. In Fig. 4, the convergence of the total energies errors is displayed for the two different correlation factors, in comparison with the coupled-cluster singles doubles perturbative triples [CCSD(T)]-F12 method. This demonstrates that improving the correlation factor can lead to a significant speed up of the basis set convergence. Using the 17-term factor, the CBS limit results can already be reached (within errors <1 mH) using cc-pVQZ basis sets.

Symmetry is a concept of paramount importance in the description and understanding of physical and chemical processes. According to Noether’s theorem, there is a direct connection between conserved quantities of a system and its inherent symmetries. Thus, identifying them allows a deeper insight into the physical mechanisms of studied systems. Moreover, the usage of symmetries in electronic structure calculations enables a much more efficient formulation of the problem at hand. The Hamiltonian formulated in a basis respecting these symmetries has a block-diagonal structure, with zero overlap between states belonging to different “good” quantum numbers. This greatly reduces the necessary computational effort to solve these problems and allows much larger systems to be studied.

There are several symmetries that are commonly used in electronic structure calculations due to the above-mentioned benefits and their ease of implementation. Our FCIQMC code NECI is no exception in this regard.

As mentioned in Sec. I, FCIQMC is usually formulated in a complete basis of Slater determinants (SDs). SDs are eigenfunctions of the total $Ŝ$_z operator, and consequently, if the studied Hamiltonian, Ĥ, is spin-independent (no applied magnetic field and spin–orbit interaction), it commutes with $Ŝ$_z, [Ĥ, $Ŝ$_z] = 0. The conservation of the m_s eigenvalue in a FCIQMC calculation thus follows quite naturally: the initial chosen m_s sector, determined by the starting SD used, will never be left by the random excitation generation process sketched in Sec. II. No terms in the spin-conserving Hamiltonian will ever cause any state in the simulation to have a different m_s value than the initial one. As a consequence, the sampled wave function will always be an eigenfunction of $Ŝ$_z with a chosen m_s, determined at the start of a calculation.

NECI is also capable of utilizing Abelian point group symmetries, with D_2h being the “largest” spatial group (similar to other quantum chemistry packages, e.g., Molcas⁸⁸ and Molpro^89,90), momentum conservation (due to translational invariance) in the Hubbard model, and uniform electron gas calculations and preservation of the m_l eigenvalues of the orbital angular momentum operator $L^z$ (the underlying molecular orbitals have to be constructed as an eigenfunction of $L^z$⁠). This is implemented via a symmetry-conserving excitation generation step and is explained in more detail in Subsection 1 a of the  Appendix.

One important symmetry of spin-preserving, nonrelativistic Hamiltonians is the global SU(2) spin-rotation symmetry. However, despite the theoretical benefits, the total SU(2) spin symmetry is not as widely used as other symmetries, like translational or point group symmetries, due to their usually impractical and complicated implementation.

There are two kinds of implementations of the total spin conservation in our FCIQMC code NECI. One approximate is based on Half-Projected Hartree–Fock (HPHF) functions.^44,91–94 Their rationale relies on the fact that for an even number of electrons, every spin state $S$ contains degenerate eigenfunctions with m_s = 0. Using time-reversal symmetry arguments, a HPHF function can be constructed as

where $Di¯$ indicates the spin-flipped version of $Di$⁠. Depending on the sign of the open-shell coupled determinants, $Hi$ are eigenfunctions of $S^2$ with odd (−) or even (+) eigenvalue S. The use of HPHF is restricted to systems with an even number of electrons and can only target the lowest even- and odd-S state. Thus, it cannot differentiate between, e.g., a singlet S = 0 and a quintet S = 2 state.

Our full implementation of the total spin conservation is based on the graphical unitary group approach (GUGA). It relies on the observation that the spin-free excitation operators Ê_ij in the spin-free formulation of the electronic Hamiltonian,

have the same commutation relations,

as the generators of the unitary group U(n). This connection allows the usage of the Gel’fand–Tsetlin (GT) basis,^95–97 which is irreducible and invariant under the action of the operators Ê_ij, in electronic structure calculations. The GT basis is a general basis for any irrep of U(n), but Paldus^98–100 realized that only a special subset is relevant for the electronic problem (34) due to the Pauli exclusion principle. Based on Paldus’s work, Shavitt¹⁰¹ further developed an even more compact representation by introducing the graphical extension of the UGA. This leads to the most efficient encoding of a spin-adapted GT basis state (CSF) in the form of a step-vector $d$⁠. This step-vector representation has the same storage cost of two bits per spatial orbital as Slater determinants. The entries of this step-vector encode the change of the total number of electrons ΔN_i and the change of the total spin ΔS_i of subsequent spatial orbitals i. This is summarized in Table I. All possible CSFs for a chosen number of spatial orbitals N, number of electrons n, and total spin S are then given by all step-vectors $d=d1,d2,…,dN$ fulfilling the restrictions,

The last restriction in Eq. (36) corresponds to the fact that the (intermediate) total spin must never be less than 0.

The most important finding of Paldus and Shavitt^102,103 was that the Hamiltonian matrix elements—more specifically the coupling coefficients between two CSFs, e.g., $m′Êijm$—can be entirely formulated within the framework of the GUGA, without any reference and thus necessity to transform to a Slater determinant-based formulation. Although CSFs can be expressed as a linear combination of SDs, the complexity of this transformation scales exponential with the number of open-shell orbitals of a specific CSF.¹⁰⁴ Thus, it is prohibitively hard to rely on such a transformation, and for already more than ≈15 electrons, a formulation without any reference to SDs is much more preferable.

Furthermore, Shavitt and Paldus^102,103 were able to find a very efficient formulation of the coupling coefficients as a product of terms via the graphical extension of the UGA. Matrix elements between two given CSFs only depend on the shape of the loop enclosed by their graphical representation, as depicted in Fig. 5. The coupling coefficient of the one-body operator Ê_ij is given by

where the product terms depend on the step-values of the two CSFs, $dk′$ and d_k, the difference in the current spin ΔS_k (with the restriction $Sk′$ − S_k = ±1/2), and the intermediate spin S_k of $m$ at orbital k. Q_k in Eq. (37) depends on the shape of the loop formed by $m$ and $m′$ at level k and was tabulated in, e.g., Ref. 102. Additionally, the two CSFs, $m$ and $m′$⁠, must coincide outside the range (i, j) for Eq. (37) to be non-zero.

The compact representation of spin-adapted basis functions in the form of step-vectors and the product form of the coupling coefficients (37) allow for a very efficient implementation in our stochastic FCIQMC code NECI. As mentioned in Sec. II, the excitation generation step is at the heart of any FCIQMC code.

The main difference to a SD-based implementation of FCIQMC, apart from the more involved matrix element calculation (37), is the higher connectivity within a CSF basis. For a given excitation operator Ê_ij, with spatial orbital indices (i, j), there is usually more than one possible excited CSF $m′$ when applied to $m$⁠, $Êijm=∑kckmk′$⁠. All valid spin-recouplings within the excitation range (i, j) can have a non-zero coupling coefficient as well. This fact is usually the prohibiting factor in spin-adapted approaches. However, there is a quite virtuous combination of the concepts of FCIQMC and the GUGA formalism, as one only needs to pick one possible excitation from $m$ to $m′$ in the excitation generation step of FCIQMC (see Sec. II).

We resolved this issue, by randomly choosing one possible valid branch in the graphical representation, depicted in Fig. 5, for randomly chosen spatial orbital indices i, j(, k, l). Additionally, we weight the random moves according to the expected magnitude of the coupling coefficients^8,105 to ensure p_gen(m′|m) ∝ |H_m′m|. This approach avoids the possible exponential scaling as a function of the open-shell orbitals of the connected states within a CSF-based approach.

However, this comes with the price of reduced generation probabilities and consequently a lower imaginary time step, as mentioned in Sec. II. Combined with an additional effort of calculating these random choices in the excitation generation and the on-the-fly matrix element computation, the GUGA-FCIQMC implementation has a worse scaling with the number of spatial orbitals N compared to a Slater determinant-based implementation.⁸ 

However, the benefits of using a spin-adapted basis are a reduced Hilbert space size, elimination of spin-contamination in the sampled wave function and, most importantly, the spin-adapted FCIQMC implementation via the GUGA allows targeting specific spin states, which are otherwise not attainable with a SD-based implementation, as discussed in Ref. 8.

The unique specification of a target spin allows resolving near degenerate spin states, and consequently, the numerical results can be interpreted more clearly. This enables more insight in the intricate interplay of nearly degenerate spin states and their effect on the chemical and physical properties of matter.

The GUGA-FCIQMC method has been benchmarked¹⁰⁵ by applying it to a linear chain of L equidistant hydrogen atoms¹⁰⁶ recently studied to test a variety of quantum chemical methods,¹⁰⁷ which shall serve as an example here. Using a minimal STO-6G basis, there is only one orbital per H atom and the system resembles a one-dimensional Hubbard model^41,108–110 with long-range interaction. Studying a system of hydrogen atoms removes complexities such as core electrons or relativistic effects and thus is a convenient benchmark system for quantum chemical methods.

For large equidistant separation of the H atoms, a localized basis, obtained with the default Boys-localization in Molpro’s LOCALI routine, with singly occupied orbitals centered at each hydrogen is more appropriate than a HF basis. Thus, this is an optimally difficult benchmark system of the GUGA-FCIQMC method, since the complexity of a spin-adapted basis depends on the number of open-shell orbitals, which is maximal for this system. Particularly, targeting the low-spin eigenstates of such highly open-shell systems poses a difficult challenge within a spin-adapted formulation. This situation is depicted schematically in Fig. 6.

We studied this system to show that we are able to treat systems with up to 30 open-shell orbitals with our stochastic implementation of the GUGA approach.¹⁰⁵ We calculated the S = 0, 1, and 2 (only S = 0 for L = 30) energy per atom up to L = 30 H atoms in a minimal STO-6G basis at the stretched r = 3.6 a₀ geometry¹⁰⁷ and compared it with the DMRG^{107,111–114} reference results. The results are shown in Table II, where we see excellent agreement within chemical accuracy with the reference results.

An important fact is the order of the orbitals though. Similar to the DMRG method, it is most beneficial to order the orbitals according to their overlap, since the number of possible spin recouplings depends on the number of open shell orbitals in the excitation range. If we make a poor choice in the ordering of orbitals, excitations between physically adjacent and thus strongly overlapping orbitals are accompanied by numerous possible spin-recouplings in the excitation range, if stored far apart in the list of orbitals. This behavior is thoroughly discussed in Ref. 115.

When applying for access to large computing clusters, it is often necessary to demonstrate that the software being used (in this case NECI) is capable of using the hardware efficiently. Ideally, the speed-up relative to using some base number of compute cores should grow perfectly linearly with the number of cores. In 2014, Booth et al.⁹⁴ presented an example with 500 × 10⁶ walkers in which no deviation from a linear speed-up is notable when compared using 512 cores to using 32, and even at 2048 cores, a speed-up by a factor of 57.5 was reported, which is 90% of the ideal speed-up factor of 64. In that work, the largest number of cores explored was 2048. By comparing the performance for a calculation with 100 × 10⁶ walkers and 500 × 10⁶ walkers, the same figure showed that the speed-up became closer to the ideal speed-up when the number of walkers was increased, suggesting that when using even more walkers, the efficiency comes even closer to 100% of the ideal speed-up factor.

Since 90% of the ideal speed-up factor was achieved in 2014 with only 500 × 10⁶ walkers on 2048 cores, and large compute clusters nowadays tend to have tens of thousands of cores available, we report scaling data for a much larger number of walkers on up to 24 800 cores in Table III. The calculations were done using the integrals in FCIDUMP format for the (54e, 54o) active space first described in Ref. 116 for the FeMoco molecule, and the output files are provided in the supplementary material.¹¹⁷ 

The scaling analysis presented in Table III was done with 32 × 10⁹ walkers on each of the two replicas used for the RDM sampling. Calculations at 32 × 10⁹ walkers are expensive, so we only completed enough iterations to determine an accurate estimate of the average runtime per iteration for the scaling analysis, and not enough iterations to accurately estimate the energy.

One may ask whether or not the scaling observed in Table III was performed for a reasonable number of walkers for this active space. To answer this question, we compare in Table IV the best (non-extrapolated) DMRG and sHCI-PT2 energies in the literature¹¹⁸ to energies obtained with i-FCIQMC at only 8 × 10⁹ walkers/replica and find that the i-FCIQMC-RDM and i-FCIQMC-PT2 energies are closer together than the sHCI-VAR and sHCI-PT2 energies, indicating that the i-FCIQMC energies are closer to the true FCI limit where the difference between variational and PT2 energies should vanish. The DMRG result lies about half-way between the two i-FCIQMC results, but fairly well below the lower of the sHCI results (a forthcoming publication specifically about the FeMoco system is planned in which more details will be presented, but the purpose of this paper is to give an overview of the NECI code).

Furthermore, comparing the time per iteration between 8 × 10⁹ and 32 × 10⁹ walkers shows that a high parallel efficiency is also achieved at a lower walker number. The determinants in NECI are stored using a hash table, making i-FCIQMC linearly scaling in the walker number,⁹⁴ so the ideal time per iteration with 32 × 10⁹ walkers at 19 960 cores according to the result for 8 × 10⁹ walkers at 16 000 cores would be 23.4 s, which is only marginally smaller than the reported 23.5 s. Note, however, that this is the relative efficiency between large scale calculations, which demonstrates performance gain from extending parallelization at large scales, not from parallelization over the entire range of scales, which is addressed to some extent by the chromium dimer example below.

In the case of the chromium dimer (cc-pVDZ, 28 electrons correlated in 76 spatial orbitals) considered in Fig. 7, the average time per iteration per walker ranges from 3.18 × 10⁻⁹ s at 640 cores to 2.51 × 10⁻¹⁰ s at 10 240 cores and 1.53 × 10⁻¹⁰ s at 20 480 cores, corresponding to a parallel speed-up of 82.1% from 10 240 to 20 480 cores and an overall speed-up of 65.2% over the full range. The deviation from ideal scaling almost exclusively stems from the communication of the spawns, at lower walker numbers, the communicative overhead is more significant, reducing the parallel efficiency compared to the FeMoco example. Nevertheless, a very high yield can be obtained from scaling up the number of cores, even for already large scales.

The parallel efficiency of NECI is made possible by treating static load imbalance. NECI contains a load-balancing feature,²⁸ which is enabled by default and periodically re-assigns some determinants to other processors in order to maintain a constant number of walkers per processor. As shown in Fig. 7 for the given benchmarks, no significant load imbalance occurs up to (including) 20 480 cores.^119,120 The initialization of a simulation does not feature the same speed-up due to I/O operations and initial communication such as trial wave function setup and core space generation. However, since it does not play a significant role for extended calculations, we consider only the time spent in the actual iterations.

The ongoing development of NECI is focused on an efficiently scaling solver for the CI-problem. It is not desirable to reimplement functionality that is already available in the existing quantum chemistry codes. Since the CI-problem is defined by the electronic integrals and subsequent methods depend on the results of the CI-step, namely, the reduced density matrices, it is easily possible to replace a CI-solver for the existing quantum chemistry code with NECI.

To use NECI, only an input file and a FCIDUMP file,¹²² which is the widely understood file format for the electronic integrals, are required. After running NECI, the stochastically sampled reduced density matrices are available as an input for further calculations in other codes. It is possible to link NECI as a library and call it directly or to run it as an external process and do the communication with explicit copying of files. The first alternative will be referred to as embedded; the second is the decoupled form.

Due to the stochastic nature of the Monte Carlo algorithm, it is not yet possible to use NECI as a black box CI-solver for larger systems. In this case, it is recommended to use the decoupled form for a better manual control of the convergence. Another advantage of the decoupled form is the combination of NECI with different quantum chemical algorithms or implementations that do not benefit from massive parallelization as much as NECI. This way it is possible to switch from serial or single node execution to multiple nodes in the CI-step. So far, NECI has been coupled with Molpro,^89,123 Molcas 8,⁸⁸ OpenMolcas,¹²⁴ PySCF,¹²⁵ and VASP.¹²⁶ 

The stochastic multi-configurational self-consistent field (MCSCF) procedure emerges from the combination of conventional MCSCF methodologies with FCIQMC as the CI-eigensolver. Stochastic-MCSCF approaches greatly enlarge the applicability of FCIQMC to strongly correlated molecular systems of practical interest in chemical science.

To date, two implementations of stochastic-MCSCF have been made available based on the interface of NECI with OpenMolcas^19,124 (and Molcas 8⁸⁸) and PySCF.^20,125 As they are both based on the complete active space (CAS) concept, they are also often referred to as stochastic-CASSCF methods.

The stochastic-CASSCF implemented in PySCF is based on a second-order CASSCF algorithm,¹²⁷ which decouples the orbital optimization problem from the active space CI problem, allowing for easy interfacing with NECI.

At each macro-iteration, a FCIQMC simulation is performed at the current point of orbital expansion, and density matrices are stochastically sampled (see Sec. VII). These are then passed back to PySCF, which updates the orbital coefficients accordingly, using either a 1-step¹²⁷ or 2-step approach.¹²⁸ 

The stochastic-CASSCF implemented in OpenMolcas is based on the quasi-second-order super-CI orbital optimization. Optimal orbitals (in the variational sense) are found by solving the super-CI secular equations in the $Super−CI$ basis, defined by the CAS wave function at the point of expansion, |0⟩, and all its possible single excitations,

The wave function is improved by mixing single excitations to the |0⟩ wave function. As the CASSCF optimization proceeds, the χ_pq coefficients decrease until they vanish, and |0⟩ will reveal the variational stationary point. Third-order density matrix elements of the exact super-CI approach are avoided by utilizing an effective one-electron Hamiltonian, as discussed in more detail in Ref. 19.

A flow chart of stochastic-CASSCF describing the various steps of the CASSCF wave function optimization is given in Fig. 8.

The stochastic-CASSCF approach has successfully been applied to a number of challenging chemical problems. The accuracy of the method has been demonstrated on simple test cases, such as benzene and naphthalene²⁰ and more complex molecular systems, namely, coronene,²⁰ free-base porphyrin, and Mg–porphyrin.¹⁹ More recently, the method has also been applied to understand the mechanism stabilizing intermediate spin states in Fe(II)–porphyrin,^39,40 the study of a $[Fe(III)2S2(SCH3)2]2−$ iron–sulfur model system in its oxidized form,¹¹⁵ and new superexchange paths in corner-sharing cuprates.¹²⁹ 

To date, only state specific stochastic-CASSCF optimizations have been reported. However, state-average stochastic-CASSCF optimizations are a straightforward extension that can be reached by taking the advantage of the NECI capability to optimize excited state wave functions, as discussed in Sec. VIII. The stochastic-CASSCF method can also be coupled to the adaptive shift approach discussed in Sec. V with a great enhancement in performance.

With NECI, we present a state-of-the-art FCIQMC program capable of running a large variety of versions of the FCIQMC algorithm. This includes the semi-stochastic FCIQMC feature, energy estimation using trial wave functions, the stochastic sampling of the reduced density matrices, and excited state calculations. Further features of NECI’s FCIQMC implementation discussed are the real-time FCIQMC method and the adaptive shift method, as well as a spin-adapted formulation of the algorithm and support for the transcorrelated Hamiltonians. We demonstrated the scalability of the program to up to 24 800 cores, showing that the code can run efficiently on large-scale machines.

Finally, we highlighted the interoperability of NECI with other quantum chemistry software, in particular OpenMolcas and PySCF, which can be used to run stochastic-CASSCF calculations.

Example FCIQMC output files for excited state calculations (output_file_excited_state_be2_b1g.txt and stats_file_excited_state_be2_b1g.txt) and real-time calculations including the resulting spectrum (output_file_real_time_be2_b1g.txt and fft_spectrum_be2_b1g.txt) for the examples presented in Sec. III are available in the supplementary material. Furthermore, it contains the output files for scaling (output_file_scaling_with_*_cores.txt and output_file_energy_with_8b_walkers.txt) and load imbalance analysis (output_file_load_imbalance_n*.txt) and also exemplary output and integral files for a similarity transformed FCIQMC calculation of the neon atom in a cc-pVDZ basis (tcdump_Ne_st_pVDZ.h5 and FCIDUMP_Ne_st_pVDZ integral files and output_file_Ne_st_pVDZ.txt and stats_file_Ne_st_pVDZ.txt files). All output files contain the corresponding FCIQMC input.

The early development of NECI was supported by the EPSRC under Grant Nos. EP/J003867/1 and EP/I014624/1.

We would like to thank Olle Gunnarsson, David Tew, Daniel Kats, Aron Cohen, and Vamshi Katakuri for insightful discussions.

The high performance benchmarks discussed in Sec. XIII were run on the MPCDF (Max Planck Computing and Data Facility) system Cobra. PJ was supported by the Marsden Fund of New Zealand (Contract No. MAU1604), from government funding managed by the Royal Society Te Apārangi. G.H.B. gratefully acknowledges the funding received from the European Research Council under Grant Agreement No. 759063.

The data that support the findings of this study are available within the article (and its supplementary material). The NECI program can be obtained at https://github.com/ghb24/NECI_STABLE, and the development version can be obtained from the corresponding author upon reasonable request.

In this appendix, we will consider in some detail the process of (random) excitation generation in FCIQMC—a crucial yet rather flexible aspect of the algorithm. We will consider some general aspects, such as implementation of Abelian symmetries in the excitation process and non-uniform excitation generation, as is often desirable in quantum chemical Hamiltonians. There are other classes of systems (such as Hubbard models, transcorrelated Hamiltonians, and spin models, such as Heisenberg systems) for which more specialized considerations are necessary for efficient excitation generation, but we will not consider them here.

The first general point about excitation generation (by which we mean starting from a given determinant $Di$⁠, we randomly pick either one or two electrons and a corresponding number of holes to substitute them with, to create a second determinant $Dj$⁠) is that if |H_ij| > 0, then the probability [p_gen(j|i)] to select $Dj$ and $Di$ must also be greater than 0. Furthermore, p_gen(j|i) must be computable, and in general, the effort to do so will depend on the algorithm chosen to execute the excitation process.

Let us discuss in more detail the process of stochastic excitation generation, and its impact on p_gen. Suppose we are simulating a system of n electrons in 2N spin orbitals {ϕ₁, …, ϕ_2N}. A given determinant $Di$ can be defined by its occupation number representation, $I=n1,…,n2N$⁠, which is a binary string such that n_i = 1, if orbital i is occupied (“an electron in $Di$”), and n_i = 0, if it is unoccupied (“a hole in $Di$”). Each orbital carries a spin quantum number σ(ϕ_i) and may also carry a symmetry label, Γ(ϕ_i). These are both discrete symmetries, with σ = ±1/2, and Γ = Γ₁,…, Γ_G, where G is the number of irreducible representations available in the point-group of the system under consideration. We will only consider Abelian groups so that the product of symmetry labels uniquely specifies another symmetry label. This simplifies the task of selecting excitations, although it does not necessarily exploit the full symmetry of the problem.

Now we wish to perform a stochastic excitation generation, which we will initially consider without the use of any symmetry/spin information. For example, we can select a pair of electrons, i, j (with i < j) in $Di$⁠, at random, and a pair of holes a, b (with a < b) and perform the transition ij → ab. The corresponding matrix element is

We will denote the electron-pair simply as ij and the hole-pair as ab.

For this simple procedure, it is clear that the probability to choose $Dj$ from $Di$ is simply

from which it follows that p_gen(j|i) ∼ (nN)⁻². This procedure does not take symmetry or spin quantum numbers into account, and it is quite possible that the corresponding Hamiltonian matrix element is zero. To ensure that we do not generate such excitations, we need to select the hole pairs so that the following two conditions are met:

These restriction greatly impact the way in which we will select i, j and a, b, and the resulting generation probability.

One way to impose symmetries in excitation generation while keeping track of the generation probabilities is via the notion of conditional probabilities. For example, rather than drawing (ij) and (ab) independently, with probability p(ab, ij) = p(ab)p(ij), one can instead draw (ab) given that one has already drawn (ij); the probability for this process is given by

where p(ij) is the probability to select (ij) in the first place. If (ij) has a particular characteristic that confers a physical (e.g., symmetry-related) constraint on (ab), this can be implemented at the stage in which we select (ab):(ab) need only be selected from among those hole-pairs for which the constraint is satisfied. For example, if the electrons (ij) have opposite spins, then the holes (ab) must also have opposite spins. The smaller number of possibilities in choosing the ab pair then leads to a larger p(ab|ij) compared to p(ab), which can be thought of as a renormalization of the latter probability to take into account the constraint.

The concept of conditional probabilities can be further extended so that the pair (ij) itself is made to satisfy a particular condition. Suppose we introduce a set of conditions ${C1,C2,…}$ such that the union of all such conditions is exhaustive. It is possible to draw conditional probabilities with respect to such conditions. For example,

and then, one can write

with

Here, $p(C1)$⁠, the probability to select the same-spin excitations, can be chosen arbitrarily, which then fixes $p(C2)$ according to the above.

The advantage of this formulation is that we can skew the selection of electron-pairs, for example, toward opposite spin excitations if that proves advantageous, and to be able to compute the resulting probabilities. Furthermore, we can write

which allows us to select a pair of electrons satisfying condition $C1$ and, subsequently, draw a pair of holes, given that one has selected an electron-pair with the same spin (which implies that the hole-pair must be chosen to have the same spin as the electron-pair).

Let us now consider how to generate the hole-pairs in a non-uniform manner, to reflect the fact that, in ab initio Hamiltonians, the matrix elements vary strongly in magnitude. Since the spawning probability is proportional to the ratio |H_ij|/p_gen(j|i), it is clearly desirable to generate excitations that make this ratio as uniform as possible, ideally with p_gen(j|i) ∝ |H_ij|. In this way, one would ensure a relatively uniform probability of successful spawning, which ideally would be close to one, implying a low rejection rate. Keeping the discussion focused on double excitations (the generalization to single excitations being straightforward), the question that arises is: how best can one select ij and ab such that p_gen(j|i) ∝ |H_ij| to a good approximation and p_gen remains exactly computable without excessive cost. We will see that there is a compromise to be made. One can ensure precise proportionality between p_gen(j|i) and |H_ij|, but only at prohibitive cost. Alternatively, one might be able to select ij and ab to affect the transition $Di→Dj$ based on computationally inexpensive heuristics to provide approximate proportionality, which will nevertheless allow for a large overall improvement in efficiency.

To ensure exact proportionality between p_gen(j|i) and |H_ij|, it is necessary to enumerate all electron-pairs and hole-pairs, which are possible from $Di$⁠, and to construct the cumulative probability function (CPF) from which the desired distribution can be straightforwardly sampled. The (unnormalized) CPF is