Automatic derivation of many-body theories based on general Fermi vacua

The formalism of second quantization, first developed in the context of the quantum field theory,¹ plays a fundamental role in many-body perturbation theory,² coupled cluster (CC) methods,^3–5 and Green’s function approaches^6–12 for molecules and solids. Two important tools for the manipulation and simplification of expressions involving second-quantized operators are Wick’s theorem¹³ and diagrammatic methods.¹⁴ These are traditionally formulated using particle/hole quasiparticle operators defined with respect to a single determinant Fermi vacuum (reference state). However, even though these tools can simplify the derivation of equations for many-body theories, this process is prone to human error and may be quite lengthy or even impossible to complete in a reasonable amount of time.

In the past three decades, computer-aided derivation of many-body equations has played an ever-increasing role in quantum chemistry.¹⁵ Starting with the pioneering works of Janssen and Schaefer¹⁶ and Li and Paldus,^17,18 the early 2000s saw the rapid development of automatic derivation and implementation tools based on algebraic,^19–26 diagrammatic,^27–30 and string or determinant-based methods.^31–35 Recently, automatic derivation has been extended to many new directions, including arbitrary-order response and derivatives,^36,37 systems with coupled fermionic and bosonic degrees of freedom,³⁸ more general vacua (Hartree–Fock–Bogoliubov, Bardeen–Cooper–Schrieffer, or antisymmetrized geminal power states),^39,40 nuclear structure theory,^41–43 and the manipulation of quantum circuits.⁴⁴ Several works have also investigated the problem of optimizing tensor operations (factorization, global optimization of the contraction order, the identification of reusable intermediates, and the identification of common factors).^45–48 

In the case of multireference (MR) many-body methods,^49–51 automatic derivation may be crucial to achieving a numerical implementation, as in many cases the underlying equations are too complicated for manual derivation and implementation. Examples of automated implementations of multireference theories include equation-of-motion multireference coupled cluster (MRCC),²⁵ arbitrary-order Mukherjee’s state-specific MRCC,⁵² internally contracted MRCC,^53–57 and analytic energy gradients and derivative couplings of multireference perturbation theories.^58,59

An important tool for the derivation of equations of internally contracted multireference theories^60–66 is the generalized normal ordering (GNO) approach proposed by Mukherjee,⁶⁷ which extends the definition of normal-ordered operators and Wick’s theorem to general correlated Fermi vacua and later was formulated using spinfree operators by Mukherjee and Kutzelnigg.⁶⁸ The GNO formalism provides the theoretical framework for extending Fock-space or “many-body equations” methods^69–72 to the multireference case. Later works have corroborated^73,74 the original proof given in Ref. 67 and examined the issue of spin adaptation.^75,76 The GNO approach is foundational to several methods, including canonical transformation theory,⁶⁵ explicitly correlated basis set incompleteness corrections,^77,78 various equation-of-motion multireference coupled cluster methods,^25,79–82 and similarity renormalization group approaches.^83–85 

This paper describes an efficient implementation of Wick’s theorem for general Fermi vacua that combines a compact representation of operator contractions with a recursive backtracking algorithm. Our representation of operators is related and generalizes diagrammatic approaches used by others.^27,54,56 The algorithm described in this work is implemented in the Wick’s theorem and diagrammatic code Wick&d (which we pronounce “wicked”), an open-source package developed by us and available from GitHub.⁸⁶ Wick&d offers similar functionality to other software designed for the derivation of multireference theories, including Neuscamman’s SQA package,^87,88 Valeev’s SeQuant2,⁸⁹ Kong’s automatic program generator APG,⁹⁰ Köhn’s GeCCo program,^56,91 Shiozaki’s Smith,^30,92 and the ADG program.^41–43 However, some of the features unique to Wick&d include the ability to generate both projective and Fock-space (many-body) equations and the support for an arbitrary number of orbital subspaces.

This article is organized as follows. In Sec. II, we summarize the generalized normal ordering formalism and the corresponding Wick’s theorem. Section III describes a general strategy to implement the generalized Wick’s theorem. Example applications of Wick&d to single-reference and multireference methods are reported in Sec. IV. In Sec. V, we conclude with a discussion of the main features, limitations, and future extension of Wick&d.

In this section, we summarize the main results of the generalized normal ordering formalism and the corresponding Wick’s theorem.^67,68,73 We follow the convention of writing products of second-quantized operators using the compact notation

where upper (lower) indices correspond to creation (annihilation) operators and are read from left to right (right to left).

Consider an arbitrary N-electron reference state Ψ₀. Mukherjee⁶⁷ defines a normal-ordered operator product ${ârs⋯pq⋯}$ to satisfy the condition

This definition can be applied recursively—starting from single substitution operators ${âsp}$—to express normal-ordered operators in terms of bare operators and reduced density matrices of the reference state $γrs⋯pq⋯=⟨Ψ0|ârs⋯pq⋯|Ψ0⟩$⁠.⁷⁴ A product of second-quantized operators $q̂1q̂2⋯$⁠, with $q̂i∈{âp†}∪{âp}$⁠, can be expressed as a sum of normal-ordered terms using a generalization of Wick’s theorem,

When compared to Wick’s theorem for a Slater determinant reference (e.g., see Refs. 14 and 93), the generalized form contains two new aspects. First, pairwise contractions yield elements of the one-particle $γ1$ or one-hole (η₁) density matrices,

Second, in addition to pairwise contractions, new multi-legged contractions appear. A 2k-leg contraction (k ≥ 2) involves k creation and k annihilation operators and corresponds to elements of the k-body density cumulant (λ_k) of the reference Ψ₀.^67,94,95 For example, the following four-leg contraction evaluates to an element of the two-body density cumulant $λrspq$⁠,

Multi-leg contractions are antisymmetric with respect to permutations of the operators involved, for example,

Two important simplifications apply to complete active space (CAS) references. First, since cumulants are zero when one or more indices correspond to empty or fully occupied spin–orbitals, multi-leg contractions only connect operators labeled by active indices. Second, the one-particle density matrix is block diagonal, with the occupied (unoccupied) block equal to the identity (zero) matrix.

A second Wick theorem helps express a product of normal-ordered operators ${Â}{B̂}⋯{Ẑ}$ as a single normal-ordered product ${ÂB̂⋯Ẑ}$ plus a sum of contractions,

Since the starting operators are normal ordered, Eq. (8) excludes contractions that only involve second-quantized operators within a normal-ordered group. This second form of Wick’s theorem plays a central role in the derivation of expressions of internally contracted multireference theories.

In this section, we describe a general procedure to evaluate Eq. (8) when the orbital space is partitioned into an arbitrary number of subspaces. We begin by defining three types of orbital subspaces and then proceed to describe a canonical representation of operators, a directed hypergraph (diagrammatic) representation of Wick contractions, and a backtracking approach for generating all unique Wick contractions.

To describe the structure of the general reference state Ψ₀, we partition the set of orthonormal spin–orbitals $S={ψ1,ψ2,…}$ into s disjoint sets $Sk$ (orbital subspaces) such that

The structure of the reference Ψ₀ is defined by constraints on the occupation of the orbitals in each subspace. We consider three types of orbital subspaces:

1. Occupied. All spin–orbitals in this subspace are occupied by one electron.

2. General. These spin–orbitals are partially occupied in the reference, and consequently, the corresponding density matrices and cumulants are nontrivial.

3. Unoccupied. All spin–orbitals in this subspace are empty.

The restrictions on the one-body density matrix, hole density matrix, and cumulants for each of these subspaces are reported in Table I. In Wick&d, the reference state is specified by the number of orbital subspaces and their type. This information is provided by the user, and it is fully customizable. For example, a mean-field reference wave function (a single Slater determinant) is specified by partitioning $S$ into occupied (⁠$O$⁠, occupied) and virtual (⁠$V$⁠, unoccupied) orbitals. A CAS reference is instead specified by splitting $S$ into core (⁠$C$⁠, occupied), active (⁠$A$⁠, general), and virtual (⁠$V$⁠, unoccupied) sets.

The automatic enumeration of Wick contractions benefits from expressing all normal-ordered operators in a canonical form. Then, an operator can be identified uniquely by the number of second-quantized operators that create and annihilate particles in each orbital subspace. To express this canonical form in a compact way, we first introduce a convenient notation for products of second-quantized operators. We write a product of $nk+$ second-quantized creation operators that act on subspace $Sk$ as

where $Pk=(p1,…,pnk+)$ is a multi-index. Similarly, we define a product of $nk−$ annihilation operators in subspace $Sk$ as

where $Qk$ is the multi-index $Qk=(q1,…,qnk−)$⁠. This notation allows us to define a product of second-quantized operators in canonical order as

where groups of creation (annihilation) operators are ordered according to increasing (decreasing) subspace index.

Using the notation of Eq. (12), we write a normal-ordered operator $Ω̂$ in a canonical form as

where $ωP1P2⋯Q1Q2⋯$ is a tensor antisymmetric with respect to separate permutations of upper and lower indices, while the numerical prefactor accounts for equivalent contributions. The term in Eq. (13) corresponds to a vertex labeled by a label (“ω”) and an operator matrix N = [n⁺n⁻], where $n+=[n1+,…,ns+]$ and $n−=[n1−,…,ns−]$ are column vectors that define the number of creation and annihilation operators in each orbital subspace, respectively. In this article, we represent such an object with a matrix and a label,

Next, we provide some examples to illustrate how this canonical representation of operators works. For a CAS reference, $S$ is partitioned into core, active, and virtual sets, $S=C∪A∪V$⁠. Then, a single excitation operator that promotes an electron from a core to an active orbital is represented by a 3 × 2 matrix,

where the corresponding second-quantized operators and entries in the matrix representation are indicated with the same color. Note that we denote this operator with the orbital subspace labels $(AC)$ corresponding to the sequence of second-quantized operators $(âu†âm)$⁠. The two-electron operator corresponding to the replacement ψ_vψ_u → ψ_mψ_e (with $ψu,ψv∈A$⁠, $ψm∈C$⁠, and $ψe∈V$⁠) is represented as follows:

where $vmeuv=〈me‖uv〉$ is an antisymmetrized two-electron integral in physicist notation. The prefactor 1/2 accounts for the equivalent contributions of terms in which the indices u and v are exchanged.

The same notation may be extended to define a product of operators. In this case, we just arrange the matrices according to the order of the operators and assign a numerical prefactor to the entire expression. For example, the product $V̂CVAA12T̂AC2$ is represented as an ordered list of operator matrices and labels multiplied by the scalar factor 1/2,

Having defined a canonical representation for normal-ordered operators, we proceed to define a canonical representation of Wick contractions [Eq. (8)]. As shown in Fig. 1(a), a Wick contraction expressed in an algebraic form may be represented as a diagram, that is, a graph in which vertices correspond to normal-ordered operators and edges represent contractions. In this interpretation, contractions (edges) can connect an arbitrary number of operators and they encode the type of second-quantized operator connected (creation or annihilation); therefore, one can establish a correspondence between diagrams and directed hypergraphs [Fig. 1(b)]. One way to represent a directed hypergraph is via an incidence matrix [Fig. 1(c)] that encodes how the edges (contractions) connect to the vertices.

When discussing Wick contractions, we distinguish between elementary (connected) and composite (disconnected) contractions. Elementary contractions are individual 2k-leg contractions involving two or more second-quantized operators and will be denoted as $Ci$⁠. To represent contractions, we use a notation similar to the one used for normal-ordered operators. For each operator, we represent an elementary contraction with a matrix C = [c⁺ c⁻], where $c+=[c1+,…,cs+]$ and $c−=[c1−,…,cs−]$ are column vectors that define the number of creation and annihilation operators contracted in each orbital subspace, respectively. Then, an elementary contraction can be represented as a list of contraction matrices,

For example, consider the single-reference partitioning $S=O∪V$ and the product

The single contraction between the leftmost creation operator $âk†$ (red) and the rightmost annihilation operator $âi$ (blue) is represented in the following way:

This representation does not specify the operators contracted within each group and may be used to designate all contractions that yield the same algebraic term. For example, in addition to the contraction shown in Eq. (20), one should also consider the following contraction obtained by connecting $âj†$ instead of $âi†$⁠,

However, as shown above, since the tensor $tabij$ is antisymmetric, after swapping the indices i and j and rearranging, this second contraction is identical to the one in Eq. (20). To account for these two equivalent contractions, when translating the incidence matrix in Eq. (20) to its algebraic form, we multiply it by a combinatorial factor of 2. This is an example of a Wick contraction that connects equivalent second-quantized operators, defined as operators of the same type (creation/annihilation) acting on the same orbital subspace. The representation of composite contractions used in Wick&d exploits this equivalence to minimize the number of terms generated and to facilitate the recognition of identical terms.

A term resulting from Wick’s theorem corresponds to a combination of elementary contractions, which we refer to as a composite contraction. A composite contraction is represented by stacking rows of contraction matrices in the hypergraph incidence matrix. Since an elementary contraction may appear more than once in a composite contraction, and the order of these is immaterial, the latter may also be represented with a multiset [e.g., ${C1,C1,C2}$]. The following example illustrates the representation of a pair of two-leg contractions:

The top matrix row is a contraction that connects a pair of virtual annihilation/creation operators (orange/blue), while the middle matrix row is a contraction involving occupied creation/annihilation operator pair (red/green). A more complex example is reported in Fig. 1, where in a CAS setting an elementary contraction of four operators yields a two-body cumulant (indicated in red).

Note that the representation adopted here may be redundant if the only valid elementary contractions are those between operators acting on the same subspace. However, there are cases when such an assumption is too restrictive. For example, one way to generate spin-integrated equations is to split the orbital subspaces into alpha and beta spin sets. Then, a term involving a mixed spin case cumulant would require expressing contractions of both alpha $(Aα)$ and beta $(Aβ)$ active orbitals, as shown in the following example:

To identify all valid operator contractions efficiently, Wick&d uses a backtracking algorithm. Subsections III D 1–III D 4 describe the steps of this algorithm.

In the first step, the code identifies all the elementary contractions between the groups of second-quantized operators. We define these as all contractions of pairs of operators for occupied/unoccupied orbital subspaces $$, $$ and the 2k-leg contractions for general orbital subspaces [see Eq. (6)]. In defining elementary contractions, we consider only one out of the potentially many equivalent contractions that can be obtained by permuting operators of the same type (creation/annihilation) that act on the same orbital subspace. For example, consider a CAS reference with orbitals partitioned as $S=C∪A∪V$⁠, and the product $V̂CCAAT̂AACC$⁠,

If we ignore the orbital labels, this term may be written as ${âC†âC†âAâA}{âA†âA†âCâC}$⁠. There are three elementary contractions for this case. A contraction of a pair of core creation-annihilation operators, represented by the following edge in the incidence matrix:

a contraction between a pair of active annihilation-creation operators,

and a four-leg contraction among all the operators that act on the active orbitals,

We employ a backtracking algorithm to identify all multisets that correspond to allowed combinations of elementary contractions. This algorithm essentially visits all the relevant branches of a tree whose leaves represent all possible multisets of elementary contractions. The backtracking algorithm used in Wick&d is both efficient and flexible since it applies to an arbitrary number of orbital subspaces. As an example, Fig. 2 shows the steps taken by the backtracking algorithm to find all composite contractions that arise from the term in Eq. (24). We start from the fully uncontracted term, represented by the empty multiset {}. Then, in step (1), we test the solution obtained by adding the first elementary contraction $(C1)$⁠, and because this is a valid contraction, we add it to the current solution, obtaining the composite contraction ${C1}$⁠. In step (2), the current contraction is combined with the elementary contraction $C1$ again, leading to ${C1,C1}$⁠, which is validated and added to the list of composite contractions. At this point, it is no longer possible to add more $C1$ contractions, and we proceed by adding the next contraction in the list, $C2$⁠. This contraction can be added up to two times, yielding the composite contractions ${C1,C1,C2}$ and ${C1,C1,C2,C2}$ [steps (3) and (4)]. At the end of step (4), all second-quantized operators are contracted, and the algorithm backtracks to ${C1,C1}$⁠. In step (5), we test and add the contraction $C3$ and then backtrack up to the contraction ${C1}$⁠. The algorithm proceeds for six more steps, at which point all valid contractions have been enumerated.

Up to this stage, we have selected only one of the possible permutations of contractions of equivalent set of operators (same type and orbitals space). Therefore, the backtracking algorithm produces a list of valid contractions that are equivalent to diagrams with distinct connectivity. Nevertheless, it is still possible to generate isomorphic diagrams that yield equivalent algebraic terms, and these are dealt with in the next step.

An important task in automatic derivation of many-body equations is expressing equivalent terms into a canonical form so that they can be collected. In Wick&d, we employ early canonicalization, whereby contractions are canonicalized before converting them to algebraic expressions. Currently, Wick&d implements an exhaustive (combinatorial) canonicalization procedure of composite contractions. Consider, for example, the commutator $[V̂AAAA,T̂AAAA]$⁠, where

If we consider the product $V̂AAAAT̂AAAA$⁠, one of the contributions to the fully contracted terms is (omitting the rows corresponding to the $C$ and $V$ spaces)

The equivalent contraction for the product $T̂AAAAV̂AAAA$ is given by

These two terms yield identical algebraic expressions, and this can be shown by expressing both of them in a canonical form. In manipulating the incidence matrix, we are allowed to reorder the contractions (rows) and operators (columns) as long as the resulting expression is equivalent to the original contraction. In the present example, starting from the contraction of $T̂AAAAV̂AAAA$⁠, we may permute the two groups of operators (by swapping columns)

and then interchange the order of contractions (by swapping the top two rows) to obtain the same incidence matrix for the corresponding term arising from $V̂AAAAT̂AAAA$⁠. Once expressed in a canonical form, these two contributions cancel, giving a zero contribution for the term $⟨Ψ0|[V̂AAAA,T̂AAAA]|Ψ0⟩$⁠.

To automate the identification of equivalent terms, we define an ordering of the directed hypergraph incidence matrices, and we define the canonical form as the minimal element out of all possible incidence matrices that represent the same term. The details of this procedure are discussed in  Appendix A.

After all composite contractions are generated and equivalent contributions are combined, each term is converted to an algebraic expression. Consider, for example, the following contraction:

To convert it to an algebraic expression, we first assign distinct indices to the operators in the order in which they appear in the second-quantized operators,

Next, we reorder this term so that second-quantized operators that appear in the same contraction are adjacent, keeping track of sign factors that arise from permutations,

After elimination of the Kronecker delta factors, we arrive at the expression

which is brought into a canonical form by rearranging the second-quantized operators and relabeling the indices. In the last step, this contribution is multiplied by a combinatorial factor (2) that keeps into account the following identical contribution that differs by a permutation of two equivalent operators:

The resulting canonical expression for the contraction ${C1,C1}$ is

The equation for the combinatorial factor of a general contraction is reported in  Appendix B. The results of Wick’s theorem lead to a sum of normal-ordered operators that may be subsequently processed or implemented as tensor contractions.

The backtracking algorithm discussed here is one the many possible ways to implement the generalized Wick’s theorem [Eq. (8)] and was chosen due to the simplicity of the recursive algorithm that generates composite contractions. We also explored directly enumerating all unique diagrams, as done by Kállay and co-workers.^27,28 However, the resulting algorithm is too complex and challenging to generalize beyond the case of a pair of operators. An alternative to the approach taken here is to build contractions of two operators at a time and apply this procedure recursively to a general product of operators on the fly. Enumerating all contractions for a pair of operators is easier, and one may be able to even directly produce only valid contractions, leading to potential speedups. However, this approach may generate an unnecessarily large number of intermediate terms at each stage of the recursion. For example, if one is interested only in the fully contracted terms for a product of operators, a pairwise recursive approach might need to store a large number of partially contracted terms.

The algebra of second-quantized operators and the backtracking algorithm to generate Wick contractions are implemented in Wick&d, an open-source C++ library exposed as a Python module via the pybind11 library.⁹⁶ Wick&d’s C++ library defines various classes used to represent diagrams (using the hypergraph incidence matrix) and algebraic terms that represent equations in terms of explicit orbital indices. The Wick&d repository includes several Jupyter-notebook tutorials on the use of the application programming interface (API) and various example applications (including the ones described in this paper). Currently, Wick&d supports the derivation of spin–orbital or spin-integrated expressions. Spin adaptation of these equations is an extension planned for future releases of the code.

In this section, we showcase two applications of Wick&d to the derivation of single-reference and multireference many-body theories.

In our first example, we derive expressions for the coupled cluster (CC) theory residuals $(rab⋯ij⋯)$⁠,

where $H̄=exp(−T̂)Ĥexp(T̂)$ is the similarity-transformed Hamiltonian. Here, the operator $T̂=T̂1+T̂2+⋯+T̂n$ is a sum of particle–hole excitation operators up to order n, with a generic operator $T̂k$ defined as

Note that our index notation for the cluster amplitudes $tab⋯ij⋯$ corresponds to the traditional notation^14,93 with upper/lower indices swapped. We obtain the CC amplitude expressions by first evaluating the many-body operator $H̄$ as follows:

and then by extracting the tensor components $H̄ab⋯ij⋯$ corresponding to the particle–hole excitation operators ${âij⋯ab⋯}$⁠. From Wick’s theorem, consequently, the kth order CC residuals $⟨Φ|{âab⋯ij⋯}H̄|Φ⟩$ are related to the tensor elements $H̄ab⋯ij⋯$ by separate antisymmetrization of the upper and lower indices,

where $Apq⋯$ is the antisymmetrizer, defined as the sum over the k! permutations of the indices p, q, … divided by 1/k!. Note that Eq. (38) applies only in the case of single-reference theories, and in the multireference case the residual is a more complicated expression of the tensor $H̄ab⋯ij⋯$⁠, density matrices, and cumulants.

The numbers of unique terms (diagrams) in CC with excitation level n ranging from 2 (CCSD) to 8 (CCSDTQPH78) computed with Wick&d are reported in Table II. At each truncation level, we confirmed that Wick&d yields the same number of diagrams reported by Kállay and Surján in their study of arbitrary order CC methods.²⁷ Together with the number of diagrams, we also report the runtime of Wick&d, showing that derivation of even CC equations with up to octuple excitations can be performed in less than a minute.

As a final test, we numerically validated the CCSD and coupled cluster with full triple excitations (CCSDT) residual equations derived by Wick&d by implementing them in a pilot Python code. We used integrals computed with Psi4⁹⁷ and generated code that calls the tensor contraction function einsum implemented in the numpy library⁹⁸ and verified that this implementation matches reference energies.

The second example uses Wick&d to implement the second-order driven similarity renormalization multireference perturbation theory (DSRG-MRPT2).⁹⁹ For the sake of brevity, we focus on the essential aspects of DSRG-MRPT2 relevant to the derivation of the corresponding equations and direct the interested reader to consult a recent review for further details.⁸⁵ In the unitary DSRG formalism, the bare Hamiltonian $(Ĥ)$ is diagonalized by a continuous unitary transformation,

where $H̄(s)$ is the transformed Hamiltonian and $Â(s)=T̂(s)−T̂†(s)$ is an anti-Hermitian operator that depends on a time-like parameter s defined in the range [0, ∞). The DSRG energy is given by the expectation value of $H̄(s)$ with respect to a complete active space self-consistent field (CASSCF) reference state Ψ₀, $E(s)=⟨Ψ0|H̄(s)|Ψ0⟩$⁠. The operator $Â(s)$ is obtained by solving a set of Fock-space many-body equations of the form $[H̄(s)]N=R̂(s)$⁠, where the subscript “N” indicates the non-diagonal part of $H̄(s)$⁠. The tensor component of the source operator $R̂(s)$ takes the form

and $rij⋯ab⋯(s)=[rab⋯ij⋯(s)]*$⁠. In Eq. (40), $Δab⋯ij⋯$ is a generalized Møller–Plesset denominator defined in terms of the diagonal components of the Fock matrix $(ϵp=fpp)$⁠,

DSRG-MRPT2 uses a diagonal Fock partitioning, whereby the Hamiltonian normal-ordered with respect to Ψ₀ is split according to $Ĥ=Ĥ(0)+ξĤ(1)$⁠. The zeroth-order Hamiltonian $Ĥ(0)$ is the sum of the reference energy (E₀) and a diagonal one-body operator $[F̂(0)]$⁠,

The first-order amplitudes that enter the operator $Â(1)(s)$ are obtained by solving the following linear equation (omitting the variable “s”):

while the second-order energy is given by

Listing 1 shows how to use Wick&d to evaluate Eqs. (43) and (44) in terms of molecular integrals and the first-order amplitude equations. In this listing, we define the $C$⁠, $A$⁠, and $V$ orbital subspaces (lines 4–6), specify the form of the operators $Ĥ(0)$⁠, $Ĥ(1)$ (lines 9–12) and the anti-Hermitian operator $Â(1)$ (lines 15–17), form the operator $H̄(1)$ (line 20), define the second-order energy (line 23), and finally apply Wick’s theorem to evaluate $H̄(1)$ and the second-order energy (lines 27 and 28). The resulting second-order energy expression contains 226 contractions. This number may be reduced to 24 contractions by introducing the following intermediate $H̃(1)=Ĥ(1)+R̂$⁠, which allows us to rewrite E⁽²⁾ as⁹⁹ 

where $H̃(1)$ is a general two-body operator,

The energy expression derived by Wick&d using Eq. (45) is reproduced as follows:

In this expression, we use Einstein notation and assign the labels to the orbital subspaces in the following way: $i,j∈C$⁠, $u,v,w,x,y,z∈A$⁠, and $a,b∈V$⁠. We numerically validated both versions of the DSRG-MRPT2 equations [Eqs. (44) and (45)] by implementing them in a pilot Python code using an interface to Forte.¹⁰⁰ This implementation also used Wick&d to derive and test equations for the first-order amplitude equations.

We have presented a computational strategy to evaluate the generalized version of Wick’s theorem^67,68 in the case of an arbitrary number of orbital subspaces. Our approach represents Wick contractions using directed hypergraphs incidence matrices. To rapidly evaluate all Wick terms, it enumerates all elementary contractions among the operators to be contracted and uses a backtracking algorithm to generate combinations of elementary contractions (composite contractions). This algorithm forms the basis for the automated derivation of many-body equations implemented in the open-source software library Wick&d. We have illustrated the utility of this code by deriving the single-reference coupled cluster equations through octuple excitations and the energy expressions for the second-order multireference driven similarity renormalization group perturbation theory (DSRG-MRPT2). The correctness of these equations was verified via pilot numerical implementations of CCSD, CCSDT, and DSRG-MRPT2.

In its current form, Wick&d is capable of deriving many-body equations that can be used for the purpose of exploring new theories and developing pilot implementations. In the case of multireference theories, Wick&d is particularly useful for the derivation of equations of internally contracted methods. Although, in our examples, we have emphasized the formalism for a general reference state, our framework is equally applicable to derive equations for state-averaged theories of multiple electronic states^101,102 and, more generally, ensembles of states. Given its ability to describe an arbitrary number of fermionic spaces, another unexplored application of Wick&d would be the derivation of many-body equations for systems involving electrons and fermionic nuclei.^103,104

In the future, we plan to expand the capabilities of Wick&d in several directions. To generate efficient implementations, it would be desirable to develop post-processing modules to spin adapt the equations generated by Wick&d and to identify reusable tensor intermediates. It would also be interesting to extend Wick’s theorem kernel to more general references. Examples include the ability to consider number-symmetry-broken vacua that require the inclusion of anomalous contractions $$e.g., among operators of the same type like $$⁠. Similarly, it would be desirable to expand Wick&d to mixed fermionic/bosonic fields to enable the derivation of many-body equations for strongly interacting electron–phonon¹⁰⁵ and electron–photon systems.^106,107

The author would like to thank Jonathon Misiewicz and Ilias Magoulas for providing feedback on a draft of this work. Thanks also go to Matthias Hanauer, Andreas Köhn, Liguo Kong, Chenyang Li, and Luther Blissett for insightful discussions in the early stages of conceptualizing the design of Wick&d. This work was supported by a grant from the U.S. Department of Energy under Award No. DE-SC0016004 and a Camille Dreyfus Teacher-Scholar Award (No. TC-18-045).

The authors have no conflicts to disclose.

The data that support the findings of this study are openly available in Wick&d, at https://github.com/fevangelista/wicked.

In this appendix, we define an ordering on the set of all possible equivalent incidence matrices representing a Wick contraction of K operators and L elementary contractions. We then discuss how this ordering is used to define a canonical form of the directed hypergraph representation. An incidence matrix $W$ may be represented by a table of matrices and labels as follows:

where N_i and ω_i are the operator matrix and label corresponding to the ith operator, respectively, while C_ij is the contraction matrix for the ith operator and jth elementary contraction.

Now consider an equivalent incidence matrix $W′$ obtained by permuting the rows and columns of the incidence matrix shown in Eq. (A1), which we write as

For two given incidence matrices $W$ and $W′$⁠, we test for the condition $W<W′$ using lexicographic ordering. Specifically, we first compare the operator labels and matrices and identify the first place i where the elements (N_i, ω_i) and $(Ni′,ωi′)$ differ. If $ωi<ωi′$⁠, then $W<W′$⁠. If $ωi=ωi′$⁠, then we test for $Ni<Ni′$⁠, by comparing the entries of the two columns of N_i and $Ni′$ lexicographically. If $ωi>ωi′$⁠, then $W′<W$⁠. If the operator matrices and labels are identical, we then proceed to compare the contraction matrices C_ij and $Cij′$⁠, one row at a time. The first place where two contraction matrices differ, we compare C_ij and $Cij′$ lexicographically. If $Cij<Cij′$⁠, then $W<W′$⁠, otherwise $W′<W$⁠.

Once an ordering of the incidence matrices is defined, we consider the subset of valid incidence matrices out of the K!L! possible ones obtained by permuting the position of the K operators and L elementary contractions. The element of this subset $W$ that is minimal (i.e., $W<W′$ for each $W′$⁠) is the Wick contraction canonical form.

In this appendix, we report the equations for the combinatorial factor associated with the incidence matrix representation of a Wick contraction. To illustrate how the combinatorial factor is derived, we use a simple example. Consider the following three-leg (red) and two-leg (blue) contractions of a group of n second-quantized operators as follows:

We assume that the operators $q̂1⋯q̂n$ are of the same type (creation or annihilation) and act on the same orbital subspace. The number of permutations of the contraction legs is given by the number of ways one can assign the three legs of the red contraction (n choose 3) times the number of ways one can assign the two legs of the blue contraction (n − 3 choose 2) to the remaining n − 3 uncontracted operators. It is easy to see that this quantity corresponds to a multinomial factor

where the multinomial factor is defined as

The total combinatorial factor for a Wick contraction is the product of combinatorial factors for each type of operator (creation/annihilation) and orbital subspace. In the case of a Wick contraction involving two or more equivalent contractions, an additional numerical factor must be included to avoid double counting contractions that are indistinguishable. The following example shows the case of a composite contraction with three equivalent elementary contractions,

which reduces the total combinatorial factor by 1/3!.

In the most general case, we consider a contraction involving K operators characterized by operator matrices $Ni=[ni+ni−]$ with $ni+=(n1i+,…,nsi+)$ and $ni−=(n1i−,…,nsi−)$⁠, where i = 1, …, K labels operators and s is the number of orbital subspaces. The composite contraction connecting these operators is defined by the contraction matrices $Cij=[cij+cij−]$⁠, where $cij+=(c1,ij+,…,cs,ij+)$ and $cij−=(c1,ij−,…,cs,ij−)$ are the number of contractions of operator i with the elementary contraction j = 1, …, L in each subspace. Then, the combinatorial factor for a contraction takes the form