Theory and implementation of a novel stochastic approach to coupled cluster

The goal of quantum chemistry is to provide accurate and cost-effective methodologies for the solution of the molecular electronic Schrödinger equation. One needs to be able not only to reproduce experimentally measurable observables but also to understand the microscopic origin of these measurements and eventually predict and guide experiments.

Stochastic approaches to the solution of the Schrödinger equation provide an appealing alternative to deterministic strategies, and a number of Monte Carlo (MC) sampling methods have been continuously developed since the early days of quantum chemistry.^1,2 At the cost of introducing statistical uncertainty in the results, quantum Monte Carlo (QMC) offers a low-scaling, parallelizable route to high-accuracy results. Despite the favorable scaling and scalability, QMC suffers from two well-known problems. The statistical errorbar can be decreased but at a very slow rate with increasing length of the simulation, that is, a larger number of random samples. Furthermore, for fermionic systems of interest in molecular electronic structure theory, the nodal structure of the wavefunction needs to be fixed a priori to avoid collapse onto the bosonic ground state,^3,4 introducing an uncontrolled approximation that, thus, far cannot be efficiently relaxed to exactness. Despite this, their low polynomial scaling allows large-scale applications, particularly within condensed matter systems where they can provide accurate results while allowing extrapolation to remove finite-size errors.^1,2,5

Within the deterministic realm, the two above-mentioned problems do not appear. No statistical uncertainity riddles the results, and the methods are all formulated in the appropriate Fock space, guaranteeing that the solution, while approximate, is properly antisymmetrized. The coupled cluster (CC) wavefunction ansatz arguably provides the most effective framework for accurate simulations of single-reference molecular systems. The CC model provides an exponential parametrization of the molecular electronic wavefunction and enjoys a number of favorable properties. It provides a systematic route toward the exact, full configuration interaction (FCI) solution, while maintaining size-extensivity and size-consistency of results at any truncation level. Despite the exponential, nonlinear parametrization of the wavefunction, the computational cost of CC theory scales as a polynomial of system size, albeit with potentially high values for the exponents. Scaling and scalability are thus much less favorable than with stochastic approaches: a number of approximations have to be introduced,^6–27 and many technical challenges need to be surmounted.^28–36 In particular, while many high-performance implementations of coupled cluster with single and double (CCSD) substitutions and CCSD with perturbative triples correction [CCSD(T)] are nowadays available, the large gain in efficiency seen in local theories has yet to be reproduced for higher truncation levels in the CC hierarchy.

With these considerations in mind, efforts in the past decade have been directed at combining the best of both worlds into the formulation and implementation of Fock-space QMC methods. The full configuration interaction quantum Monte Carlo (FCIQMC) was the first such method to be presented, and the FCI secular problem is solved as the dynamics of a population of signed particles.^37–40 This results in an exponentially scaling algorithm but with a dramatically reduced prefactor. Building on this, some of us have further developed a similar projector MC algorithm to solve the unlinked and linked CC equations.^41–45 These CCMC algorithms, implemented in the HANDE-QMC software package,⁴⁶ are fully general with respect to the excitation level, allowing one to perform arbitrary order CC simulations with a sparse representation of the wavefunction.

We recently showed that neither FCIQMC nor CCMC rigorously fulfills size-extensivity for noninteracting systems.⁴⁷ Both algorithms perform imaginary-time propagation of unlinked many-body equations,⁴⁸ which results in the unnecessary sampling of zero-on-average terms. This unnecessary work negatively impacts the memory and central processing unit (CPU) costs of the simulation and is particularly severe for CCMC, as it quickly precludes scaling to larger systems and/or higher orders of CC theory. To remedy this situation, we put forth a MC algorithm that performs the imaginary-time propagation governed by the linked CC equations. These are evaluated by the random sampling of the connected terms in the similarity-transformed Hamiltonian, conveniently represented as a diagrammatic expansion. The diagCCMC algorithm restores size-extensivity, and our preliminary tests have shown how localization can be readily exploited without further assumptions.⁴⁹ 

We should note that this is not the only avenue toward leveraging the benefits of MC sampling within the CC approach. Deustua et al. have shown how deterministic iterative CC solvers can be seeded with amplitudes from partially converged Fock-space QMC results. Combined with moment expansion corrections,^50–52 this approach is a powerful technique, enabling access to higher levels of CC theory at a reduced cost. However, the high computational scaling of the QMC methods used to determine important higher-level amplitudes will eventually dominate the overall computational cost of this approach, and thus, our work provides a complementary solution.

In this work, we will first describe in detail the theoretical framework in which our diagCCMC algorithm rests. Section II summarizes CC theory with particular emphasis on its diagrammatic formulation. In Sec. III, we present a derivation of the imaginary-time update step and its usage in the CCMC and diagCCMC algorithms. We will then discuss the structure of the implemented algorithm and highlight differences and similarities to a deterministic implementation of CC theory.

We will use the tensor notation for second quantization.^53–55 We denote the elementary, anticommuting fermion creation and annihilation operators as:

A k-electron excitation operator with respect to the physical vacuum (|vac⟩) is the product of k creation and k annihilation operators. In tensor notation,

such excitation operators are particle number-conserving. Explicitly, the one- and two-electron substitutions are

The Born–Oppenheimer molecular electronic Hamiltonian is then expressed as

where the integrals are given in an orthonormal basis of one-electron spin–orbitals,

For single-reference theories, it is more convenient to work in terms of the Fermi, rather than the physical, vacuum state. Our Fermi vacuum will be a single-determinant reference function |D₀⟩ for an N-electron system with 2M spin–orbitals.

Occupied one-particle states in |D₀⟩, i₁, i₂, …, i_N, will be referred to as hole states, whereas virtual one-particle states, a_N+1, a_N+2, …, will be referred to as particle states. A normal-ordered k-electron substitution operator will be denoted as

Using Wick’s theorem,^55,56 these operators can be rewritten as a finite sum of subsets of permutations of elementary operators times contractions. The latter are elements of k-electron reduced density matrices (k-RDMs),

Thus, the normal-ordered one- and two-electron substitutions are⁵⁷ 

Imposing normal ordering on the molecular Hamiltonian, we obtain

where the energy of the reference determinant, the one-body Fock operator, and the two-body fluctuation potential appear,

We will use the symbol $τ^k$ for pure excitation operators or excitors. These are k-electron substitutions between hole and particle states in the reference determinant and, thus, particle number- and charge-conserving.

Since the k-RDMs in (7) are zero whenever any of the indices, upper or lower, refers to a particle state, the excitors are automatically normal-ordered,

Here, we have introduced the multi-index $k=a1a2…aki1i2…ik$ to compactly represent the k-electron substitution affected by the excitor, which is, up to a phase, a k-excited determinant,

A summary of our notation can be found in Table I.

The coupled cluster wavefunction is parametrized as an exponential transformation of a reference single-determinant wavefunction |D₀⟩,

where the cluster operator $T^$ is given as a sum of second-quantized excitation operators,

with the mth order cluster operators expressed as sums of excitors weighted by the corresponding cluster amplitudes,

Note that in the tensor notation adopted, upper and lower indices of the cluster amplitudes appear reversed with respect to other conventions.

The CC correlation energy is the right eigenvalue of the Schrödinger equation for the normal-ordered Hamiltonian in (9),

This equation is solved by performing a similarity transformation of the Hamiltonian,

and then projecting onto the excitation manifold {|D_j⟩},

Equation (18b) defines the CC residual ω_k(t), which is zero at a solution of the nonlinear linked equations.⁵⁸ Whereas (16) and (18a) can be proved to be identical,⁵⁹ the linked formulation in the latter is size-extensive order-by-order and term-by-term. For notational convenience, we have dropped the subscript N for the Hamiltonian. The similarity-transformed Hamiltonian $H¯$ can be expanded into a Baker–Campbell–Hausdorff (BCH) commutator series,

For the molecular Hamiltonian in Eq. (4), at most two-body operators are involved. Hence, regardless of the truncation level in the cluster operator $T^$⁠, the expansion truncates at the fourfold nested commutator,

showing that only finitely many terms are included in Eq. (18b). Despite the fact that $H¯$ is no longer Hermitian, the linked formulation is still more advantageous than the unlinked formulation.^59–61 

Since all excitors are normal-ordered and commuting, Wick’s theorem^56,60,61 lets us reduce the Hamiltonian-excitor products to only those terms that are connected (in the diagrammatic sense). Excitors will only appear to the right of the Hamiltonian, and only terms where each excitor shares at least one index with the Hamiltonian will be nonzero,

The requirement of shared indices between the Hamiltonian and cluster coefficients enables the resulting equations to be solved via a series of tensor contractions: a process highly amenable to rapid evaluation on conventional computing architectures⁶² but non-trivial to parallelize.^28,29,35

The algebraic derivation of the linked CC equations to a general truncation level from Eq. (18b) is lengthy and error prone. A diagrammatic representation can be effectively used to generate all unique terms in the equations.^60,63,64 The normal-ordering and application of Wick’s theorem are key to these developments. The normal-ordered Hamiltonian features 13 interaction vertices, 4 coming from the Fock operator,

and 9 from the fluctuation potential,

Each of these vertices can be characterized by an integer representing their excitation level (0, ±1, ±2) and by a sign sequence encoding the pattern of open particle (+) and hole (−) lines below the interaction vertex (see Table II). Cluster operators can be classified similarly in terms of their excitation level (any integer ≥1) and their sign sequence.

For any given excitation level in the allowed manifold (up to double excitations for CCSD, triple excitations for CCSDT, and so forth), the diagrammatic generation of the corresponding CC equations proceeds via the following steps:

1. At the bottom, we draw a combination of at most four excitors.

2. At the top, we draw a Hamiltonian vertex. The valid vertices are limited by two requirements: (a) the final diagram be connected and (b) the overall excitation level of the projection manifold.

3. We pair the Hamiltonian vertex and excitor(s) sign sequences in all distinct ways to generate the sign sequences for all unique diagrams. The sign sequence encodes the diagram topology and ensuing contraction pattern.

4. We read the algebraic expression for the corresponding term in the CC equations off from the generated diagrams. The rules of interpretation associate target indices to the external (open) lines and dummy summation indices to the internal lines, Hamiltonian matrix elements to the interaction vertices, and products of amplitudes to the excitor vertices. Topological and permutational symmetries are taken into account by similar simple rules.^60,61,63

The rules for generating and interpreting diagrams as algebraic expressions are independent of the CC truncation order and can be encoded into a computer program.^61,65–70 However, a proper factorization of intermediates is essential to achieve acceptable time to solution and memory requirements.⁶⁶ 

The solution of the CC equations can be achieved by means of stochastic algorithms. This stochastic realization is, however, not unique, and multiple algorithms have been put forward in the literature.^41–43 All these different realizations are based on reformulating the time-dependent Schrödinger equation in imaginary time. The corresponding diffusion-like equation can be solved by the repeated application of an approximate propagator on a trial state. Employing a Fock-space representation circumvents the fermion sign problem without the need for fixing the nodes a priori.⁷¹ 

After performing a Wick rotation τ ← it to imaginary time, the time-dependent CC Schrödinger equation reads as^72,73

The τ-derivative on the left-hand side is (see the  Appendix)

Excitation operators are assumed to be time-independent,

and since all excitors commute, the nested commutator expansion truncates at l = 0,

The imaginary-time Schrödinger equation (24) then becomes

and upon projection onto $Dk|exp−T^(τ)=D0|τ^k†⁡exp−T^(τ)$⁠,

since by construction $⟨D0|τ^k†τ^j|D0⟩=δkj$⁠. Equation (29) is an imaginary-time ordinary differential equation (ODE) that we can solve by discretization.

The stochastic propagation of the linked CC equations is, thus, directly related to those utilized within FCIQMC,³⁷ diffusion Monte Carlo (DMC),^1,74 and the original unlinked coupled cluster Monte Carlo (CCMC) approach.^41,42 This allows us to understand limits on the time step due to the spectral range of the Hamiltonian and more directly compare computational costs with prior stochastic coupled cluster theory.

The imaginary-time ODE in Eq. (29) can be discretized in a number of ways. In principle, we would like to (a) use as large a time step as possible without losing the stability of the integrator and (b) perform the fewest possible number of evaluations of the CC vector function per time step. The usual approach in CCMC and FCIQMC is the explicit Euler method with a time step h,

where $tk[n+1]$ and $tk[n]$ are the cluster amplitudes at times τ + h and τ, respectively, and $ωk[n]$ is the CC vector function at time τ.

Alternatively, one could use an implicit Euler scheme,

where the right-hand side now depends on the CC vector function evaluated at time τ + h. We can approximate this term using the Newton–Raphson step,⁵⁹ 

where the CC Jacobian has been introduced,

and obtain the Rosenbrock–Euler method,^75,76

Under the assumption of non-singular Jacobian, we can use a Woodbury-type identity to compute the inverse,⁷⁷ 

and retaining the first term only yields the deterministic Newton–Raphson step,⁵⁹ 

Given this point of view, it is possible to relate the imaginary-time propagation to a number of standard techniques in numerical analysis. Given a time-step δτ, the generalized step,

will be equivalent to a relaxed Newton–Raphson method.

The use of the full CC Jacobian for preconditioning would be extremely expensive, and a more pragmatic route is taken in practice. The simplest choice is to approximate the Jacobian with the identity matrix, i.e., no preconditioning is applied to the iterations. A more sophisticated approach is to only retain iteration-independent terms in Eq. (33),

where the “d” and “od” stand for diagonal and off-diagonal, respectively. We can then propose two cheap preconditioners. We can either use the diagonal part of the Fock operator,⁷⁸ 

or the diagonal part of the full Hamiltonian,

The former is universally implemented in deterministic CC codes, and its effectiveness can be justified through perturbative arguments.⁵⁹ The use of the latter has not, to the best of our knowledge, been attempted before.

The derivation here presented makes explicit connection with preconditioning, already discussed by some of us in connection with FCIQMC⁷⁹ and unlinked CCMC.⁸⁰ We will discuss how preconditioning is implemented for diagCCMC in Sec. IV C.

Finally, let us point out that Jarlebring et al. showed how a specific instance of a nonlinear eigenvalue problem is equivalent to a Rosenbrock-type discretization of an associated imaginary-time ODE.⁷⁶ An adaptive time step integrator can be, thus, formulated based on convergence estimates similar to those presented in Ref. 76.

We wish to stochastically solve the linked CC equation (18b). Additionally, and at variance with the approach of Franklin et al., we wish to overcome the need for a corrected update step and the sampling of extraneous unlinked terms.⁴³ Whereas the latter have been observed to cancel out on average, they impose limitations to what system sizes are approachable before the memory cost becomes prohibitive.

In the diagCCMC algorithm,⁴⁷ we use the uncorrected update step in Eq. (30). Two novel insights allow us to achieve this goal:

• The CC wavefunction is stored in a compressed representation without invoking particles or walkers. It is comparatively easier to enforce constant unit intermediate normalization within a walker-less algorithm.

• The CC vector function appearing in the update step is an integral expressible as a terminating series expansion. Terms in this expansion can be evaluated stochastically.

The use of diagrammatic techniques automatically guarantees that only connected terms in the similarity-transformed Hamiltonian are included. The sampling will, thus, happen in “diagram space” and relies on the even selection algorithm of Scott et al.⁴⁴ 

Previous algorithms to stochastically solve the linked CC equations modified the propagation in (30) to approach the correct solution. The need for such modifications can be attributed to the use of a variable intermediate normalization,

where the additional normalization parameter N₀ is constrained by the energy equation,

and the unknown CC energy must be substituted by the shift S. At the beginning of the simulation, S = E_ref and this causes the energy estimator to converge incorrectly prior to initialization of population control. However, upon closer inspection, the wavefunction ansatz in (41) is seen to be equivalent to the conventional CC ansatz with (a) overlap with the reference set to N₀ and (b) all nonzero cluster amplitudes t_k represented by values larger than $1N0$⁠. The floating intermediate normalization can then be interpreted as an algorithmic choice to determine the granularity of representation during the calculation and achieve compression of the CC wavefunction. This choice is arbitrary and can be related back to the conventional CC ansatz. Assume then that the intermediate normalization is now a constant value ⟨D₀|CCMC⟩ = N₀, set as an input parameter to the calculation. At sufficiently small granularities, the calculation will spontaneously stabilize at a system-dependent population of walkers without the need for population control. The stochastic realization of the modified explicit Euler integration,

would then

1. compress the cluster amplitudes to the selected granularity, by stochastically rounding those amplitudes for which $|t̃k|<1$ to $sgn(t̃k)×1$ or 0,

2. evaluate the CC vector function by taking a large enough number of samples such that diagrams in $ω̃k$ of magnitude 1 are, on average, selected once,

3. adjust the time step δτ as to avoid particle blooms, which are a large spawning event, which would destabilize the calculation dynamics.

We can, however, take one further step and cast away the walker interpretation entirely. The thresholding implied in the previous algorithmic sketch can be rigorously formulated without recourse to walkers. We introduce three strictly positive calculation parameters: the representation granularity Δ, the evaluation granularity γ, and the maximum diagram contribution ϵ. The algorithm then will

1. compress the cluster amplitudes to the chosen representation granularity, by stochastic rounding amplitudes for which |t_k| < Δ to sgn(t_k) × Δ or 0,

2. evaluate the CC vector function stochastically such that diagrams with magnitude γ are selected once on average,

3. adjust the time step δτ such that the maximum diagram contribution, $δτwdiagrampdiagram$⁠, is of magnitude ϵ.

The walker and walker-less representations are entirely equivalent. The representation granularity is the inverse of the intermediate normalization constant $Δ=1N0$⁠, the condition γ = Δ defines the even selection approach,⁴⁴ and the ratio $ϵΔ$ is the maximum allowed size for a spawning event. The resultant approach to the imposition of sparsity bears some resemblance to recent fast randomized iteration approaches.^81,82

Within this approach, the total walker population is the sum of rescaled cluster coefficient absolute magnitudes and the reference $1Δ+∑i|ti|Δ$⁠. It is, thus, not needed to set the hard-to-predict total walker population as a calculation parameter: choosing to stochastically round all cluster coefficients below a certain value gives a more intuitively stable treatment between different calculations. The total walker population can vary dramatically with system size: evaluating and comparing computational cost and performance for systems of varying size can be a nontrivial challenge. Instead, we expect the walker-less picture to manifest the transferability property of cluster amplitudes:⁸³ the magnitude of the amplitudes should be relatively unchanged with system size, especially when localized orbitals are used, allowing equivalent parameters for different calculations to be easily identified.

We have found γ = 10⁻³ to be the lowest evaluation granularity giving a calculation stable enough to extract statistics from. While smaller γ values achieve more stable calculations, with 10⁻⁴ providing a reasonable compromise between the computational cost and stability. We have also continued to use conventions from the particle representation for now by setting $Δγ=1$ and $ϵγ=3$⁠.

The second essential insight enabling the diagCCMC algorithm is the stochastic evaluation of the CC vector function on the right-hand side of the uncorrected update step. At any given excitation level in the CC hierarchy, the BCH expansion of $H¯$ will truncate at the fourfold nested commutator: ω_k is expressible as a sum of a finite enumerable number of terms. We choose to represent these terms as diagrams and generate such an expansion on the fly rather than enumerating the allowed diagrams beforehand. In each main Monte Carlo cycle in the algorithm, we perform the evaluation of the integral by attempting to select n_a fully specified diagrams from its expansion. The action of the similarity-transformed Hamiltonian on the reference determinant can be written compactly as

where the amplitude $w$_l = $w$_H ∏_mt_m is a product of a one- or two-body integral from the Hamiltonian and a cluster of excitors. The multi-index l is fully specified, meaning that all hole and particle lines are explicitly labeled. The amplitude is determined by the contraction pattern randomly selected during diagram generation. Finally, since $⟨D0|τ^p†τ^q|D0⟩=δpq$⁠, the selected diagram can contribute to one and only one cluster amplitude: the one whose multi-index corresponds to the external lines in $w$_l. The rules for the deterministic enumeration of diagrams that were briefly detailed in Sec. II C are largely unmodified in our stochastic algorithm. Each step corresponds to an event occurring with an easily computed probability:

1. Sample the action of the wave operator: with probability p_sel, choose a term from the BCH expansion (20), that is, select a cluster of size N ≤ 4 and the excitation level of each constituent excitor. We use the even selection scheme of Scott et al.⁴⁴ in a walker-less representation (see Sec. V A).

2. Sample the action of the Hamiltonian: with probability p_hver, choose one of the 13 interaction vertices in Table II. This step is not independent of the former, and we use importance sampling (see Sec. V A).

3. Sample the admissible contraction patterns: with probability p_cont, choose a specific Kucharski–Bartlett sign sequence,^60,61,63 see Table III, for example.

4. Sample the index set to label internal lines. Given the number of internal hole and particle lines in the selected contraction, the probability associated with this step p_int is computed combinatorially.

5. Sample the index set to label external lines. As for the previous step, the probability p_ext is also computed combinatorially.

With this process, we are able to obtain a given diagram with probability p_diagram = p_selectp_hverp_contp_intp_ext, and in each Monte Carlo step, the diagram is sampled p_diagram × n_a times.

We need further minor modifications to the deterministic enumeration of diagrams to ensure that the CC vector function is evaluated correctly. Our algorithm singles out specific diagrams, where all lines, internal and external, are explicitly labeled. This procedure identifies a single cluster amplitude t_k to which the selected diagram will contribute without having to sum over internal lines. Permutational symmetries will, thus, have to be handled differently such that our probability distributions are properly normalized. Sums of the form $12∑i,j$ have to be replaced with $∑i>j+12δij$ to ensure that there is only a single way to select diagrams related by the following:

• the antipermutation of indices stemming from antisymmetrized interaction vertices,

• the antipermutation of hole or particle indices stemming from excitor vertices, and

• the commutation of excitors.

For the first two cases, terms with i = j would vanish when summing over equivalent indices. In the last case, the diagonal case i = j indicates additional symmetries of the resulting diagram. In our stochastic diagram enumeration, each pair of equivalent internal or external lines will not require a $12$ factor. Moreover, upon selection, a well-defined ordering of excitors is established, which removes the need for $12$ factors in diagrams where excitors of the same rank appear. These modifications to the deterministic evaluation rules ensure a unique selection of a contraction pattern.^60,61,63 The action of permutation operators for inequivalent external lines is subsumed into the permutation of hole and particle indices needed to store the result of the sampling in antisymmetrized ordering, which provides the appropriate parity factor (−1)^σ. With these considerations, a contribution to t_k is computed as

and the sampling algorithm is designed such that p_diagram = |$w$_diagram|. Ultimately, our aim is to achieve importance sampling between contributions (see Sec. V C).

While the imaginary-time propagation discussed in Sec. III will generally be used within our work, we could also make use of arbitrary preconditioners. Apart from the identity, we implemented two additional options: the diagonal of the Fock operator (diagrams 1 and 2 in Table II) and the diagonal of the full Hamiltonian (diagrams 1, 2, 4, 5, and 6 in Table II). The connected portions of these vertices do not modify excitors when applied. These preconditioners are the iteration-independent approximations to the Jacobian discussed in Sec. III B and strike a balance between the computational complexity and improvement of convergence. The diagonal Fock preconditioner is ubiquitously implemented in deterministic CC approaches.⁵⁹ Since all relevant quantities are precomputed, the values of either preconditioner can be evaluated with a cost independent of the system size, unlike the implementation of the similar approach within FCIQMC and CCMC.^79,80

The portion of the Hamiltonian used for preconditioning can then be applied via a straightforward rescaling of the original cluster amplitudes by a factor of 1 − δτ. The remainder of the Hamiltonian is applied explicitly, as in the imaginary-time propagation, before rescaling by the preconditioner.

We will not investigate the benefits of preconditioning here but wanted to observe that the diagrammatic formalism lends itself to a straightforward implementation of a range of preconditioners without introducing additional computational costs scaling with system size. This results from the use of the connected portions of all preconditioners, unlike previous stochastic approaches.^79,80

The even selection algorithm was proposed by Scott and Thom⁴⁴ to improve the sampling of the action of the CC wave operator on the reference determinant. Even selection was specifically designed to alleviate calculation instabilities due to the occurrence of large particle blooms. Sampling proceeds via the selection of clusters containing a specific number of excitors of each rank, termed a combination, separately. In this section, we summarize the adaptation of even selection in a walker-less context. We then illustrate the need for the importance sampling of $(H⁡expT)c$ and describe the strategy implemented in diagCCMC.

The original even selection algorithm defined the probability of selecting a particular set of excitors e from the combination c of size s as

a series of conditional probabilities. We re-express the selection probability as follows: