Efficient Formulation of Full Configuration Interaction Quantum Monte Carlo in a Spin Eigenbasis via the Graphical Unitary Group Approach

Efficient Formulation of Full Configuration Interaction Quantum Monte Carlo in a Spin Eigenbasis via the Graphical Unitary Group Approach

Werner Dobrautz w.dobrautz@fkf.mpg.de Max Planck Institute for Solid State Research, Heisenbergstr. 1, 70569 Stuttgart, Germany    Simon Smart simondsmart@gmail.com European Centre for Medium-Range Weather Forecasts, Shinfield Rd, Reading RG2 9AX, United Kingdom    Ali Alavi Max Planck Institute for Solid State Research, Heisenbergstr. 1, 70569 Stuttgart, Germany Dept of Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 1EW, United Kingdom a.alavi@fkf.mpg.de
July 23, 2019

We provide a spin-adapted formulation of the Full Configuration Interaction Quantum Monte Carlo (FCIQMC) algorithm, based on the Graphical Unitary Group Approach (GUGA), which enables the exploitation of SU(2) symmetry within this stochastic framework. Random excitation generation and matrix element calculation on the Shavitt graph of GUGA can be efficiently implemented via a biasing procedure on the branching diagram. The use of a spin-pure basis explicitly resolves the different spin-sectors and ensures that the stochastically sampled wavefunction is an eigenfunction of the total spin operator . The method allows for the calculation of states with low or intermediate spin in systems dominated by Hund’s first rule, which are otherwise generally inaccessible. Furthermore, in systems with small spin gaps, the new methodology enables much more rapid convergence with respect to walker number and simulation time. Some illustrative applications of the GUGA-FCIQMC method are provided: computation of the spin gap of the cobalt atom in large basis sets, achieving chemical accuracy to experiment, and the , , , spin-gaps of the stretched N molecule, an archetypal strongly correlated system.

Quantum Monte Carlo, SU(2) symmetry
02.70.Ss, 31.10.+z, 31.15.xh, 31.25.−v

I Introduction

The concept of symmetry is of paramount importance in physics and chemistry. The exploitation of the inherent symmetries and corresponding conservation laws in electronic structure calculations not only reduces the degrees of freedom by block-diagonalization of the Hamiltonian into different symmetry sectors, but also ensures the conservation of “good” quantum numbers and thus the physical correctness of calculated quantities. It also allows to target a specific many-body subspace of the problem at hand. Commonly utilized symmetries in electronic structure calculations are discrete translational and point group symmetries, angular momentum and projected spin conservation.

Due to a non-straight-forward implementation and accompanying increased computational cost, one often ignored symmetry is the global spin-rotation symmetry of spin-preserving, nonrelativistic Hamiltonians, common to many molecular systems studied. This symmetry arises from the vanishing commutator


and leads to a conservation of the total spin quantum number .

In addition to the above-mentioned Hilbert space size reduction and conservation of the total spin , solving for the eigenstates of in a simultaneous spin-eigenbasis of allows targeting distinct—even (near-)degenerate—spin eigenstates, which allows the calculation of spin gaps between states inaccessible otherwise, and facilitates a correct physical interpretation of calculations and description of chemical processes governed by the intricate interplay between them. Moreover, by working in a specific spin sector, convergence of projective techniques which rely on the repeated application of a propagator to an evolving wavefunction is greatly improved, especially where there are near spin-degeneracies in the exact spectrum.

The Full Configuration Interaction Quantum Monte Carlo (FCIQMC) approach Booth, Thom, and Alavi (2009); Cleland, Booth, and Alavi (2010) is one such methodology which can be expected to benefit from working in a spin-pure many-body basis. Formulated in Slater determinant (SD) Hilbert spaces, at the heart of the FCIQMC algorithm is excitation generation, in which from a given Slater determinant, another Slater determinant (a single or double excitation thereof) is randomly selected to be spawned on, with probability and sign determined by the corresponding Hamiltonian matrix element. Such individual determinant-to-determinant moves cannot, in general, preserve the total spin, which instead would require a collective move involving several SDs. Therefore, although the FCI wavefunction is a spin eigenvector, this global property of the wavefunction needs to emerge from the random sampling of the wavefunction, and is not guaranteed from step to step. Especially in systems in which the wavefunctions consist of determinants with many open-shell orbitals, this poses a very difficult challenge. If, instead, excitation generation between spin-pure entities could be ensured, this would immensely help in achieving convergence, especially in the aforementioned problems.

To benefit from the above mentioned advantages of a spin-eigenbasis, we present in this work the theoretical framework to efficiently formulate FCIQMC in a spin-adapted basis, via the mathematically elegant unitary group approach (UGA) and its graphical (GUGA) extension, and discuss the actual computational implementation in depth.

There are several other schemes to construct a basis of eigenfunctions, such as the Half-Projected Hartree-Fock (HPHF) functions Smeyers and Doreste-Suarez (1973); Helgaker, Jørgensen, and Olsen (2000), Rumer spin-paired spin eigenfunctions Rumer (1932); Weyl, Rumer, and Teller (1932); Simonetta, Gianinetti, and Vandoni (1968); Smart (2013); Reeves (1966), Kotani-Yamanouchi (KY) genealogical spin eigenfunctions Kotani and Amemiya (1955); Van Vleck and Sherman (1935); Pauncz (1979), Serber-type spin eigenfunctions, Serber (1934); Pauncz (1979); Salmon and Ruedenberg (1972), Löwdin spin-projected Slater determinants Löwdin (1955) and the Symmetric Group Approach Duch and Karwowski (1982); Ruedenberg (1971); Pauncz (1995)—closely related to the UGA—, which are widely used in electronic structure calculations. Some of these have partially been previously implemented in FCIQMC (HPHF, Rumer, KY and Serber)—but with severe computational limitations. Booth and Alavi et. al. (2013); Booth et al. (2011); Booth, Smart, and Alavi (2014). The GUGA approach turns out to be quite well suited to the FCIQMC algorithm, and is able to alleviate many of the problems previously encountered.

Concerning other computational approaches in electronic structure theory, there is a spin-adapted version of the Density Matrix Renormalization Group algorithm McCulloch and Gulácsi (2002); Tatsuaki (2000); Zgid and Nooijen (2008); Sharma and Chan (2012); Li and Chan (2017), a symmetry-adapted cluster (SAC) approach in the coupled cluster (CC) theory Ohtsuka et al. (2007); Nakatsuji and Hirao (1978, 1977), where is conserved due to fully spin- and symmetry-adapted cluster operators and the projected CC method Qiu et al. (2017); Tsuchimochi and Ten-no (2019); He and Cremer (2000); Tsuchimochi and Ten-no (2018), where the spin-symmetry of a broken symmetry reference state is restored by a projection, similar to the Löwdin spin-projected Slater determinants Löwdin (1955).

The use of spin-eigenfunctions in the Columbus Lischka et al. (2001, 2011, 2017), Molcas Aquilante et al. (2016) and GAMESS software package Schmidt et al. (1993); Gordon and Schmidt (2005) packages rely on the graphical unitary group approach (GUGA), where the CI method in GAMESS is based on the loop-driven GUGA implementation of Brooks and Schaefer Brooks and Schaefer (1979); Brooks et al. (1980).

Based on the GUGA introduced by Shavitt Shavitt (1977, 1978), Shepard et al. Shepard and Simons (1980); Lischka et al. (1981) made extensive use of the graphical representation of spin eigenfunctions in form of Shavitt’s distinct row table (DRT). In the multifacet graphically contracted method Shepard (2005, 2006); Gidofalvi and Shepard (2009); Öhrn et al. (2010); Shepard, Gidofalvi, and Brozell (2014a, b); Gidofalvi, Brozell, and Shepard (2014) the ground state wavefunction is formulated nonlinearly based on the DRT, conserving the total spin .

In this paper, we begin by reviewing the GUGA approach, concentrating on those aspects of the formalism that are especially relevant to the FCIQMC method, including the concept of branching diagrams in excitation generation. We then present a brief overview of the FCIQMC algorithm in the context of the GUGA method, including a discussion of optimal excitation generation and control of the time step. Next we provide application of this methodology to spin-gaps of the N atom, the N molecule and the cobalt atom, which illustrate several aspects of the GUGA formulation. In Sec. IX we conclude our findings and give an outlook to future applications and possible extensions or our implementation.

Ii The (Graphical) Unitary Group Approach

In this section we discuss the use of the Unitary Group Approach (UGA) Paldus (1974) to formulate the FCIQMC method in spin eigenfunctions. The UGA is used to construct a spin-adapted basis—also known as configuration state functions (CSFs)—, which allows to preserve the total spin quantum number in FCIQMC calculations. With the help of the Graphical Unitary Group Approach (GUGA), introduced by Shavitt Shavitt (1977), an efficient calculation of matrix elements entirely in the space of CSFs is possible, without the necessity to transform to a Slater determinant (SD) basis. The GUGA additionally allows effective excitation generation, the cornerstone of the FCIQMC method, without reference to a non spin-pure basis and the need of storage of auxiliary information.

In this work we concern ourselves exclusively with spin-preserving, nonrelativistic Hamiltonians in the Born-Oppenheimer approximation Born and Oppenheimer (1927) in a finite basis set. The basis of the unitary group approach (UGA), which goes back to Moshinsky Moshinsky (1968), is the spin-free formulation of the spin-independent, non-relativistic, electronic Hamiltonian in the Born-Oppenheimer approximation, given as


where .With the reformulation

we can define




as the singlet one- and two-body excitation operators Helgaker, Jørgensen, and Olsen (2000), which do not change and upon acting on a state, , with definite total and z-projection value of the spin. With Eqs. (3) and (4) the Hamiltonian (2) can be expressed in terms of these spin-free excitation operators as Matsen (1964)


where . An elegant and efficient method to create a spin-adapted basis and calculate the Hamiltonian matrix elements in this basis is based on the important observation that the spin-free excitation operators (3) and (4) in the non-relativistic Hamiltonian (5) obey the same commutation relations as the generators of the Unitary Group Paldus (1974, 1975, 1976), being the number of spatial orbitals. The commutator of the spin-preserving excitation operators can be calculated as


which is the same as for the basic matrix units and the generators of the unitary group .

The Unitary Group Approach (UGA) was pioneered by Moshinsky Moshinsky (1968), Paldus Paldus (1974) and Shavitt Shavitt (1977, 1978), who introduced the graphical-UGA (GUGA) for practical calculation of matrix elements. With the observation that the spin-free, nonrelativistic Hamiltonian (5) is expressed in terms of the generators of the unitary group, the use of a basis that is invariant and irreducible under the action of these generators is desirable. This approach to use dynamic symmetry to block-diagonalize the Hamiltonian is different to the case where the Hamiltonian commutes with a symmetry operator. In the UGA does not commute with the generators of , but rather is expressed in terms of them. Block diagonalization occurs, due to the use of an invariant and irreducible basis under the action of these generators. Hence, the UGA is an example of a spectrum generating algebra with dynamic symmetry Iachello (1993); Sonnad et al. (2016).

ii.1 The Gel’fand-Tsetlin Basis

The Gel’fand-Tsetlin (GT) Gel’fand and Cetlin (1950a, b); Gel’fand (1950) basis is invariant and irreducible under the action of the generators of . The group has generators, , and a total of Casimir operators, commuting with all generators of the group, and the GT basis is based on the group chain


where is Abelian and has one-dimensional irreducible representations (irreps). Each subgroup
has Casimir operators, resulting in a total of commuting operators, named Gel’fand invariants Gel’fand (1950). The simultaneous eigenfunctions of these invariants form the GT basis and are uniquely labeled by a set of integers related to the eigenvalues of the invariants. Thus, based on the branching law of Weyl Weyl (1931), a general -electron CSF can be represented by a Gel’fand pattern Gel’fand and Cetlin (1950a)


The integers in the top row (and all subsequent rows) of (8) are nonincreasing, , and the integers in the subsequent rows fulfill the condition


called the “in-between” condition Louck (1970).

The non-increasing integers of the top row of Eq. (8), , are called the highest weight or weight vector of the representation and specify the chosen irrep of ; the following rows uniquely label the states belonging to the chosen irrep.

Let be the irreducible representation of , uniquely specified by the weight vector . Any representation of a group G yields a representation of any of its subgroups , , subduced by , . of subduced by is simply reducible Paldus (2006), due to the branching law of the unitary group Weyl (1946),


where the direct sum extends over all irreps of for which the ”in-between” condition (9) holds and each irrep is contained once at most Louck (1970); Weyl, Rumer, and Teller (1932).

This fact and since is Abelian with one-dimensional irreps led Gel’fand and Tsetlin to the realization that the permissible highest weights of the subgroups in the chain (7) can be used to uniquely label the basis vectors of a general irrep space.

In CI calculations one usually employs a one-particle basis of spin-orbitals with creation and annihilation operators of electrons in spatial orbital with spin . The operators


can be associated with the generators of with the commutation relation


The partial sums over spin or orbital indices of these operators


are related to the orbital and spin generators. (Superscript denotes the pure orbital space and the pure spin space.) Since we deal with fermions we have to restrict ourselves to the totally antisymmetric representations of , denoted as . Since the molecular Hamiltonian (5) is spin independent, we can consider the proper subgroup of the direct product of the spin-free orbital space , with generators , and the pure spin space with the four generators Paldus (1974), given as


With the total antisymmetric representation of , and and as the highest weights representing the irreps of and respectively, the subduced representation of contains only those representations of for which and are mutually conjugate Paldus (2006); Matsen (1974, 1964); Moshinsky (1968).

Plainly spoken, this means the irreps of and are related in a specific manner to obtain physically plausible states satisfying the Pauli exclusion principle and antisymmetry of fermionic wavefunctions. This is another aspect of the fact that in the totally antisymmetric wavefunction, the totally antisymmetric orbital part must be combined with the totally symmetric spin part and vice versa. E.g. an antisymmetric spin function forces a symmetric spatial function ( or ), yielding an antisymmetric singlet state. On the other hand, a symmetric spin function () is combined with an antisymmetric spatial function to yield the antisymmetric triplet states.

Moreover, since the Hamiltonian (5) is spin-independent, does not contribute to the matrix element evaluation, so we only have to concern ourselves with the irreps of the orbital subgroup, following Matsen’s spin-free approach Matsen (1974, 1964).

The consequence of the mutually conjugate relationship between and irreps for electronic structure calculations is that the integers in a Gel’fand pattern (8) for are related to occupation numbers of spatial orbitals. This means they are restricted to , due to the Pauli exclusion principle. The highest weight, , indicates the chosen electronic state with the conditions


with being the total number of electrons and the number of singly occupied orbitals, is equal to twice the total spin value .

ii.2 The Paldus table

The restriction of in electronic Gel’fand patterns led Paldus Paldus (1974) to the more compact formulation by a table of integers. It is sufficient to count the appearances , and in each row of a Gel’fand pattern and store this information, denoted by and in a table, named a Paldus table.

The first column, , contains the number of doubly occupied orbitals, the second column, , the number of singly occupied and the last one, , the number of empty orbitals, as shown by the example of an state:


where the differences , with , of subsequent rows are also indicated. For each row the condition


holds, thus any two columns are sufficient to uniquely determine the state. The top row satisfies the following properties


completely specifying the chosen electronic state as an irrep of .

The total number of CSFs for a given number of orbitals , electrons and total spin is given by the Weyl-Paldus Paldus (1974); Weyl (1946) dimension formula


As it can be seen from Eq. (19), the number of possible CSFs—of course—still scales combinatorially with the number of electrons and orbitals, as seen in Fig. 1 with a comparison to the total number of possible SDs (without any symmetry restriction). The ratio of the total number of SDs and CSFs for can be estimated by Stirling’s formula (for sufficiently large and ) as


which shows orbital dependent, , decrease of the efficient Hilbert space size for a spin-adapted basis. The Paldus table also emphasizes the cumulative aspects of the coupling between electrons, with the i-th row providing information on number of electrons, (up to i-th level) and the spin, , by


As can be seen in Eq. 16, there are four permissible difference vectors (, with ) between consecutive rows of a Paldus table, which corresponds to the possible ways of coupling a spatial orbital based on the group chain (7). This information can be condensed in the four-valued step value, shown in Table 1.

All possible CSFs of a chosen irrep can then be encoded by the collection of the step values in a step-vector, where starting from the “vacuum” ’th row , an empty spatial orbital is indicated by , a “positively spin-coupled” orbital, , by , a “negatively spin-coupled”, , by and a doubly occupied spatial orbital by . To retain physically allowed states the condition applies. (As a side note: Another common notation —e.g. in Molcas—is to indicate positive spin-coupling as , negative spin-coupling by and a doubly occupied orbital by .)

0 0 -0 1 0 0
1 0 -1 0 1 -1/2
2 1 -1 1 1 -1/2
3 1 -0 0 2 0
Table 1: Four possible ways of coupling an orbital .

The step-value in Tab. 1 is given by and the collection of all into the step-vector representation is the most compact form of representing a CSF, with the same storage cost as a Slater determinant, with 2 bits per spatial orbital. One can create all basis function of a chosen irrep of by constructing all possible distinct step-vectors which lead to the same top-row of the Paldus table (18), specifying the chosen irrep with definite spin and number of electrons, with the restriction .

Iii The Graphical Unitary Group Approach (GUGA)

The graphical unitary group approach (GUGA) of Shavitt Shavitt (1977, 1981) is based on this step-vector representation and the observation that there is a lot of repetition of possible rows in the Paldus tables specifying the CSFs of a chosen irrep of . Instead of all possible Paldus tables, Shavitt suggested to just list the possible sets of distinct rows in a table, called the distinct row table (DRT). The number of possible elements of this table is given by Shavitt (1977)


with , which is drastically smaller than the total number of possible CSFs (19) or Slater determinants (without any symmetry restrictions) as seen in Fig. 1. Each row in the DRT is identified by a pair of indices , with being the level index, related to the orbital index and being the lexical row index such that if or if and .

A simple example of the DRT of a system with , and is shown in Table 2.

a b c i j
2 0 1 3 1 2 0 3 4 - - - -
2 0 0 2 2 0 0 0 5 1 0 0 0
1 1 0 2 3 0 5 0 6 0 0 1 0
1 0 1 2 4 5 0 6 7 0 0 0 1
1 0 0 1 5 0 0 0 8 4 3 0 2
0 1 0 1 6 0 8 0 0 0 0 4 3
0 0 1 1 7 8 0 0 0 0 0 0 4
0 0 0 0 8 - - - - 7 6 0 5
Table 2: Distinct row table for , and .

Relations between elements of the DRT belonging to two neighboring levels and are indicated by the so called downward, , and upward, , chaining indices, with . These indices indicate the connection to a lexical row index in a neighboring level by a step-value , where a zero entry indicates an invalid connection associated with this step-value. Given a DRT table any of the possible CSFs can be generated by connecting distinct rows linked by the chaining indices.

Figure 1: (Color online) Number of total SDs (without any symmetry restrictions), CSFs and entries of the distinct row table (DRT) for and as a function of .
Figure 2: (Color online) Graph representing the DRT of Table 2. The orange line corresponds to the CSF and the green line to in the step-vector representation.

This DRT table can be represented as a graph, see Fig. 2, where each distinct row is represented by a vertex (node) and nonzero chaining indices are indicated by an arc (directed edge). The vertices are labeled according to the lexical row index , starting at the unique head node at the top, which corresponds to the highest row . It ends at the second unique null row , which is called the tail of the graph. Vertices with the same -value of Table 2 are at the same level on this grid. The highest -value is on top and the lowest at the bottom. Vertices also have left-right order with respect to their value and vertices that share the same value are further ordered—still horizontally—with respect to their value. With the above mentioned ordering of the vertices according to their and values, the slope of each arc is in direct correspondence to the step-value , connecting two vertices. corresponds to vertical lines, and the tilt of the other arcs increases with the step-value .

Each CSFs in the chosen irrep of , is represented by a directed walk through the graph starting from the tail and ending at the head, e.g. the green and orange lines in Fig. 2 (color online), representing the states and in step-vector representation. Such a walk spans arcs (number of orbitals) and visits one node at each level . There is a direct correspondence between the Paldus table, Gel’fand patterns and directed walks on Shavitt graphs for representing all possible CSFs in a chosen irrep of .

iii.1 Evaluation of Nonvanishing Hamiltonian Matrix Elements

Given the expression of the nonrelativistic spin-free Hamiltonian in (5) a matrix element between two CSFs, and , is given by:


The matrix elements, and , provide the coupling coefficients between two given CSFs and and are the integral contributions. The coupling coefficients are independent of the orbital shape and only depend on the involved CSFs, and - Therefore, for a given set of integrals the problem of computing Hamiltonian matrix elements in the GT basis is reduced to the evaluation of these coupling coefficients. The graphical representation of CSFs has been proven a powerful tool to evaluate these coupling coefficients thanks to the formidable contribution of Paldus, Boyle, Shavitt and others Paldus and Boyle (1980); Shavitt (1978); Downward and Robb (1977).

The great strength of the graphical approach is the identification and evaluation of nonvanishing matrix elements of the excitation operators (generators) , between two GT states (CSFs), . The generators are classified according to their indices, with being diagonal weight (W) and with being raising (R) and lowering (L) operators (or generators). In contrast to Slater determinants, applied to yields a linear combination of CSFs ,


with an electron moved from spatial orbital to orbital without changing the spin of the resulting states . They are called raising (lowering) operators since the resulting will have a higher (lower) lexical order than the starting CSF .

The distance, , from to , is an important quantity and is called the range of the generator . For the one-body term in (5) Shavitt Shavitt (1977) was able to show that the walks on the graph, representing the CSFs and , must coincide outside of this range to yield a non-zero matrix element. The two vertices in the DRT graph, related to orbital and (with ) represent the points of separation of the walks and they are named loop head and loop tail. And the matrix element only depends on the shape of the loop formed by the two graphs in the range , shown in Fig. 3.

Figure 3: (Color online) Graphical representation of a matrix element as a loop shape created by two CSFs and on a Shavitt graph.

Shavitt Shavitt (1978) showed that the relations


between and must be fulfilled to yield a nonzero matrix element ( for a raising and for a lowering generator).

This allows two possible relations between the vertices at each level in terms of Paldus array quantities depending on the type of generator (R,L). For raising generators R:


where and for lowering generators L:


At each vertex of the loop in range one of the relations (26-29) must be fulfilled for the one-body matrix element to be non-zero.

Based on the graphical approach, Shavitt Shavitt (1978) showed that the matrix elements of the generators can be factorized in a product, where each term corresponds to a segment of the loop in the range and is given by


where is the value of state at level . additionally depends on the segment shape of the loop at level , determined by the type of the generator , the step values and and . The nonzero segment shapes for a raising (R) generator are shown in Fig. 4. In Table 3 the nonzero matrix elements of the one-electron operator —an over/under-bar indicates the loop head/tail—depending on the segment shape symbol, the step-values and the -value are given in terms of the auxiliary functions

Figure 4: Nonzero segment shapes of a raising generator . The numbers next to the lines indicate the step-values and . correspond to the loop tail (head) segments and to shapes inside the generator range . indicates the possible difference of and leading to nonzero matrix elements.
00 0 01 1 1 10 1 1
11 1 02 1 1 20 1 1
22 1 13 31
33 2 23 32
00 1 1 1 1
11 -1 -1
12 - - -
21 - - -
22 -1 -1
33 -1 -1 -1 -1
Table 3: Nonzero matrix elements of the one-body operator in terms of the auxiliary functions (31.

iii.2 Two-Body Matrix Elements

The matrix elements of the two-body operators are more involved than the one-body operators, especially the product of singlet excitation generators, . Similar to the one-electron operators, the GT states and must coincide outside the total range to for to be nonzero. The form of the matrix element depends on the overlap range of the two ranges


One possibility to calculate the matrix element would be to sum over all possible intermediate states, ,


but in practice this is very inefficient. For non-overlapping ranges the matrix element just reduces to the product


where must coincide with in the range and with in range . The same rules and matrix elements as for one-body operators apply in this case. An example of this is shown in the left panel of Fig. 5.
For , we define the non-overlap range


where the same restrictions and matrix elements as for one-body operators apply. In the overlap range, , different restrictions for the visited Paldus table vertices apply for the matrix element to be nonzero. This depends on the type of the two generators involved and were worked out by Shavitt Shavitt (1981). For two raising generators (RR) the following conditions apply


For two lowering generators (LL):


And for a mixed combination of raising and lowering generators (RL)


Drake and Schlesinger Drake and Schlesinger (1977), Paldus and Boyle Paldus and Boyle (1980), Payne Payne (1982) and Shavitt and Paldus Shavitt (1981) were able to derive a scheme, where the two-body matrix elements can be computed as a product of segment values similar to the one-body case (30)


where and are the overlap (32) and non-overlap (35) ranges defined above.
are the already defined single operator segment values, listed in Table 3, and are new segment values of the overlap range (their listing is omitted for brevity here, but can be found in Refs. [dccclxxx(69; 74)]. The sum over two products in corresponds to the singlet coupled intermediate states (), with a nonzero contribution if and the triplet intermediate coupling ().

Figure 5: (Color online) Two examples of nonzero two-body matrix elements. (a) shows a non-overlapping () raising-lowering (RL) example and (b) shows a overlapping () loop with two raising generators (RR).

This product formulation of the two-body matrix elements in a spin-adapted basis is the great strength of the graphical unitary group approach, which allows an efficient implementation of the GT basis in the FCIQMC algorithm. The details of the matrix element calculation in this basis are, however, tedious and will be omitted here for brevity and clarity. More details on the matrix element calculation, especially the contributions of the two-body term to diagonal and one-body matrix elements can be found in Appendix B or in Refs. [dccclxxx(69; 74)].

Iv Spin-Adapted Full Configuration Interaction Quantum Monte Carlo

The Full Configuration Interaction Quantum Monte Carlo (FCIQMC) method Booth, Thom, and Alavi (2009); Cleland, Booth, and Alavi (2010) attempts to obtain the exact solution of a quantum mechanical problem in a given single-particle basis set by an efficient sampling of a stochastic representation of the wavefunction—originally expanded in a discrete antisymmetrised basis of Slater determinants (SDs)—through the random walk of walkers, governed by imaginary-time the Schrödinger equation. For brevity of this manuscript we refer the interested reader to Refs. [dccclxxx(1; 2)] and [dccclxxx(21)] for an in-depth explanation of the FCIQMC method.

Having introduced the theoretical basis of the unitary group approach (UGA) and its graphical extension (GUGA) to permit a mathematically elegant and computationally efficient incorporation of the total spin symmetry in form of the Gel’fand-Tsetlin basis, here we will present the actual implementation of these ideas in the FCIQMC framework, termed GUGA-FCIQMC.

Fundamentally, the three necessary ingredients for an efficient spin-adapted formulation of FCIQMC are:

  1. Efficient storage of the spin-adapted basis

  2. Efficient excitation identification and matrix element computation

  3. Symmetry adapted excitation generation with manageable computational cost

The first point is guaranteed with the UGA, since storing the information content of a CSF and a SD amounts to the same memory requirement, with CSFs represented in the step-vector representation. Efficient identification of valid excitations is rather technical and explained in Appendix A and in Ref. [dccclxxx(74)]. For the present discussion we simply need to know, although it is more involved to determine if two CSFs are connected by a single application of than for SDs, it is possible to do so efficiently. Matrix element computation is based on the product structure of the one- (30) and two-body (45) matrix elements derived by Shavitt Shavitt (1978) explained above and presented in more detail in App. B and in Ref. [dccclxxx(74)]. Concerning point (iii): symmetry adaptation in FCIQMC is most efficiently implemented at the excitation generation step, by creating only symmetry-allowed excitations. For the continuous spin symmetry this is based on Shavitt’s DRT and the restriction for nonzero matrix elements in the GUGA. This, in addition to the formulation in a spin-pure GT basis, ensures that the total spin quantum number is conserved in a FCIQMC calculation.

iv.1 Excitation Generation: Singles

The concept of efficient excitation generation in the spin-adapted GT basis via the GUGA will be explained in detail by the example of single excitations. Although more complex, the same concepts apply for generation of double excitation, which are discussed below.

In contrast to excitation generation for SDs, there are now two steps involved for a CSF basis. The first, being the same as in a formulation of FCIQMC in Slater determinants, is the choice of the two spatial orbitals and , with probability . This should be done in a way to ensure the generation probability to be proportional to the Hamiltonian matrix element involved. However, here comes the first difference of a CSF-based implementation compared to a SD-based one. For Slater determinants, the choice of an electron in spin-orbital and an empty spin-orbital is sufficient to uniquely specify the excitation , and to calculate the involved matrix element . However, in a CSF basis, the choice of an occupied spatial orbital , and empty or singly occupied spatial orbital , only determines the type of excitation generator acting on an CSF basis state as well as the involved integral contributions and of the matrix element . To ensure , the occupied orbital and (partially) empty are picked in the same way as for SDs, but with an additional restriction to ensure . However, the choice of does not uniquely determine the excited CSF as there are multiple possible ones, as explained above.

As a consequence, the choice of spatial orbitals and does not determine the coupling coefficient of the matrix element . Optimally, for a given and generator , the connected CSF has to be created with a probability proportional to the coupling coefficient . By ensuring is proportional to the integral contributions and to the coupling coefficients, the total spawning probability


will be proportional to the magnitude of Hamiltonian matrix element . The efficiency of the FCIQMC algorithm depends on the ratio of the Hamiltonian matrix element between two connected states and the probability to choose the excitation , as the imaginary timestep of the simulation is adapted to faithfully account for all excitations


In a primitive implementation, is determined by the “worst-case” ratio during a simulation. A less strict approach to this problem is discussed below. By choosing nonzero and ensuring is achieved by a branching tree approach, we obtain one of the different possible walks on the Shavitt graph with nonzero loop contributions with the starting CSF .

iv.2 The Branching Tree

In the spin-adapted excitation generation, after a certain generator is picked with a probability based on the integral contributions of the Hamiltonian matrix element, the type of generator is determined, raising (R) if and lowering (L) if . One connecting single excitation is then chosen by looping from starting orbital to and stochastically choosing a valid nonzero Shavitt graph, based on the restrictions (26-29), mentioned in the GUGA section above. As an example, let us have a closer look at a chosen raising generator. As can be seen in the single segment value Table 3 there are 4 possible nonzero starting matrix elements. These starting segments are associated with a relative difference of the total spin