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

TL;DR

This paper introduces a spin-adapted FCIQMC method using GUGA, enabling efficient and accurate calculations of spin states in strongly correlated systems by exploiting SU(2) symmetry.

Contribution

The authors develop a novel spin-pure formulation of FCIQMC based on GUGA, allowing for efficient sampling of low and intermediate spin states in complex systems.

Findings

Achieved chemical accuracy for cobalt atom's spin gap in large basis sets.

Enabled rapid convergence in systems with small spin gaps.

Successfully computed spin gaps in stretched N2 molecule.

Abstract

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 $\hat{\mathbf S}^2$. 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 $^2F-^4F$ spin gap of the cobalt atom in large basis sets, achieving chemical accuracy to experiment, and the $^1\Sigma_g^+$, $^3\Sigma_g^+$, $^5\Sigma_g^+$, $^7\Sigma_g^+$ spin-gaps of the stretched N$_2$ molecule, an archetypal strongly correlated system.

Tables

Table 1. Table 1: Four possible ways of coupling an orbital i 𝑖 i .

$d_i$	$\Delta a_i$	$\Delta b_i$	$\Delta c_i$	$\Delta N_i$	$\Delta S_i$
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 2. Table 2: Distinct row table for n = 3 𝑛 3 n=3 , N = 4 𝑁 4 N=4 and S = 0 𝑆 0 S=0 .

a	b	c		i	j		$k_0$	$k_1$	$k_2$	$k_3$		$l_0$	$l_1$	$l_2$	$l_3$
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 3. Table 3: Nonzero matrix elements of the one-body operator E ^ i ​ j subscript ^ 𝐸 𝑖 𝑗 \hat{E}_{ij} in terms of the auxiliary functions ( 30 .

$d^{\prime }d$	W		$d^{\prime }d$	$\overline{R}$	$\underset{¯}{L}$		$d^{\prime }d$	$\underset{¯}{R}$	$\overline{L}$
00	0		01	1	1		10	1	1
11	1		02	1	1		20	1	1
22	1		13	$A(b,0,1)$	$A(b,2,1)$		31	$A(b,1,0)$	$A(b,0,1)$
33	2		23	$A(b,2,1)$	$A(b,0,1)$		32	$A(b,1,2)$	$A(b,2,1)$
				R				L
			$d^{\prime }d$	$\Delta b=-1$	$\Delta b=+1$			$\Delta b=-1$	$\Delta b=+1$
			00	1	1			1	1
			11	-1	$C(b,0)$			$C(b,1)$	-1
			12	-$1/\left(b+2\right)$	-			$1/\left(b+1\right)$	-
			21	-	$1/b$			-	-$1/\left(b+1\right)$
			22	$C(b,2)$	-1			-1	$C(b,1)$
			33	-1	-1			-1	-1

Table 4. Table 4: Nonzero starting segments for R ¯ ¯ 𝑅 \underline{R} with the number of electrons N i ′ = N i + 1 superscript subscript 𝑁 𝑖 ′ subscript 𝑁 𝑖 1 N_{i}^{\prime}=N_{i}+1 in all cases.

$\underset{¯}{R}$					$\overline{R}$				$\underset{¯}{L}$					$\overline{L}$
$d_i$	$d_i^{\prime }$	$\Delta S_i$	$\Delta b_i$		$d_j$	$d_j^{\prime }$	$\Delta b_{j-1}$a		$d_i$	$d_i^{\prime }$	$\Delta S_i$	$\Delta b_i$		$d_j$	$d_j^{\prime }$	$\Delta b_{j-1}$a
0	1	$+1/2$	$-1$		1	0	$-1$		1	0	$-1/2$	$+1$b		0	1	$+1$
0	2	$-1/2$	$+1$c		2	0	$+1$		2	0	$+1/2$	$-1$		0	2	$-1$
1	3	$-1/2$	$+1$b		3	1	$+1$		3	1	$+1/2$	$-1$		1	3	$-1$
2	3	$+1/2$	$-1$		3	2	$-1$		3	2	$-1/2$	$-1$c		2	3	$+1$

Table 5. Table 5: Nonzero intermediate R 𝑅 R and L 𝐿 L segments.

$d_k$	$d_k^{\prime }$	$\Delta b_{k-1}=-1$	$\Delta b_{k-1}=+1$
		$\Delta b_k$	$\Delta b_k$
0	0	$-1$	$+1$
1	1	$-1$	$+1$
1	2	$-1$a	-c
2	1	-c	$-1$
2	2	$-1$	$-1$b
3	3	$-1$	$+1$

Table 6. Table 6: Distinct types of double excitations. i < j < k < l 𝑖 𝑗 𝑘 𝑙 i<j<k<l in all cases and e i ​ j , k ​ l = e k ​ l , i ​ j subscript 𝑒 𝑖 𝑗 𝑘 𝑙 subscript 𝑒 𝑘 𝑙 𝑖 𝑗 e_{ij,kl}=e_{kl,ij} in mind.

Label	Generator order	Operator
0a	$\underset{¯}{R}\left(i\right)\to \overline{R}\left(j\right)$	${\hat{E}}_{ij}$
0b	$\underset{¯}{R}\left(i\right)\to \overline{R}\underset{¯}{R}\left(j\right)\to \overline{R}\left(k\right)$	${\hat{e}}_{ij,jk}$
0c	$W\underset{¯}{R}\left(i\right)\to \overline{R}\left(j\right)$	${\hat{e}}_{ii,ij}$
0d	$\underset{¯}{L}\left(i\right)\to \overline{L}\left(j\right)$	${\hat{E}}_{ji}$
0e	$\underset{¯}{L}\left(i\right)\to \overline{L}\underset{¯}{L}\left(j\right)\to \overline{L}\left(k\right)$	${\hat{e}}_{ji,kj}$
0f	$\underset{¯}{L}\left(i\right)\to W\overline{L}\left(j\right)$	${\hat{e}}_{ji,jj}$
0g	$\underset{¯}{L}\left(i\right)\to W\left(j\right)\to \overline{L}\left(k\right)$	${\hat{e}}_{ji,kk}$
1a	$\underset{¯}{L}\left(i\right)\to \overline{L}\underset{¯}{R}\left(j\right)\to \overline{R}\left(k\right)$	${\hat{e}}_{ji,jk}$
1b	$\underset{¯}{R}\left(i\right)\to \overline{R}\underset{¯}{L}\left(j\right)\to \overline{L}\left(k\right)$	${\hat{e}}_{ij,kj}$
1c	$\underset{¯}{R}\left(i\right)\to \underset{¯}{R}R\left(j\right)\to \overline{RR}\left(k\right)$	${\hat{e}}_{jk,ik}$
1d	$\underset{¯}{L}\left(i\right)\to L\underset{¯}{L}\left(j\right)\to \overline{LL}\left(k\right)$	${\hat{e}}_{ki,kj}$
1e	$\underset{¯}{L}\left(i\right)\to \underset{¯}{R}L\left(j\right)\to \overline{RL}\left(k\right)$	${\hat{e}}_{jk,ki}$
1f	$\underset{¯}{R}\left(i\right)\to \underset{¯}{L}R\left(j\right)\to \overline{RL}\left(k\right)$	${\hat{e}}_{kj,ik}$
1g	$\underset{¯}{RR}\left(i\right)\to R\overline{R}\left(j\right)\to \overline{R}\left(k\right)$	${\hat{e}}_{ik,ij}$
1h	$\underset{¯}{LL}\left(i\right)\to \overline{L}L\left(j\right)\to \overline{L}\left(k\right)$	${\hat{e}}_{ji,ki}$
1i	$\underset{¯}{RL}\left(i\right)\to \overline{R}L\left(j\right)\to \overline{L}\left(k\right)$	${\hat{e}}_{ij,ki}$
1j	$\underset{¯}{RL}\left(i\right)\to \overline{L}R\left(j\right)\to \overline{R}\left(k\right)$	${\hat{e}}_{ji,ik}$
2a	$\underset{¯}{RR}\left(i\right)\to \overline{RR}\left(j\right)$	${\hat{e}}_{ij,ij}$
2b	$\underset{¯}{LL}\left(i\right)\to \overline{LL}\left(j\right)$	${\hat{e}}_{ji,ji}$
2c	$\underset{¯}{RL}\left(i\right)\to \overline{RL}\left(j\right)$	${\hat{e}}_{ij,ji}$
3a	$\underset{¯}{R}\left(i\right)\to \underset{¯}{R}R\left(j\right)\to R\overline{R}\left(k\right)\to \overline{R}\left(l\right)$	${\hat{e}}_{jl,ik}/{\hat{e}}_{jk,il}$
3b	$\underset{¯}{L}\left(i\right)\to L\underset{¯}{L}\left(j\right)\to \overline{L}L\left(k\right)\to \overline{L}\left(l\right)$	${\hat{e}}_{ji,lk}/{\hat{e}}_{li,kj}$
3c0	$\underset{¯}{R}\left(i\right)\to \overline{R}\left(j\right)\to \underset{¯}{R}\left(k\right)\to \overline{R}\left(l\right)$	${\hat{e}}_{ij,kl}$
3c1	$\underset{¯}{R}\left(i\right)\to \underset{¯}{L}R\left(j\right)\to \overline{L}R\left(k\right)\to \overline{R}\left(l\right)$	${\hat{e}}_{kj,il}$
3d0	$\underset{¯}{L}\left(i\right)\to \overline{L}\left(j\right)\to \underset{¯}{L}\left(k\right)\to \overline{L}\left(l\right)$	${\hat{e}}_{ji,lk}$
3d1	$\underset{¯}{L}\left(i\right)\to \underset{¯}{R}L\left(j\right)\to \overline{R}L\left(k\right)\to \overline{L}\left(l\right)$	${\hat{e}}_{jk,il}$
3e0	$\underset{¯}{R}\left(i\right)\to \overline{R}\left(j\right)\to \underset{¯}{L}\left(k\right)\to \overline{L}\left(l\right)$	${\hat{e}}_{ij,lk}$
3e1	$\underset{¯}{R}\left(i\right)\to \underset{¯}{L}R\left(j\right)\to \overline{R}L\left(k\right)\to \overline{L}\left(l\right)$	${\hat{e}}_{lj,ik}$
3f0	$\underset{¯}{L}\left(i\right)\to \overline{L}\left(j\right)\to \underset{¯}{R}\left(k\right)\to \overline{R}\left(l\right)$	${\hat{e}}_{ji,kl}$
3f1	$\underset{¯}{L}\left(i\right)\to \underset{¯}{R}L\left(j\right)\to \overline{L}R\left(k\right)\to \overline{R}\left(l\right)$	${\hat{e}}_{jl,ki}$

Table 7. Table 7: Automatically obtained timesteps for an SD- and GUGA-based (optimized and vanilla) simulation of N 2 at equilibrium distance in a cc-pVDZ basis. Results for the “worst-case” optimization Δ ​ τ w Δ subscript 𝜏 𝑤 \Delta\tau_{w} and for the integrated histogram based optimization Δ ​ τ h Δ subscript 𝜏 ℎ \Delta\tau_{h} covering 99.99% of all excitations.

	$\Delta {\tau }_w$	$\Delta {\tau }_h$	$\Delta {\tau }_h/\Delta {\tau }_w$
SD	$5.59\cdot {10}^{-3}$	$6.20\cdot {10}^{-3}$	1.11
GUGA van.	$4.78\cdot {10}^{-5}$	$8.62\cdot {10}^{-5}$	1.80
GUGA opt.	$5.20\cdot {10}^{-5}$	$1.12\cdot {10}^{-3}$	21.50
GUGA van. / opt.	0.92	0.08
SD / GUGA opt.	107.51	5.55

Table 8. Table 8: Spin gap D o 2 − S o 4 superscript superscript 𝐷 𝑜 2 superscript superscript 𝑆 𝑜 4 {}^{2}D^{o}-{{}^{4}S^{o}} and ionization potential (IP) P 0 3 − S o 4 superscript subscript 𝑃 0 3 superscript superscript 𝑆 𝑜 4 {}^{3}P_{0}-{{}^{4}S^{o}} of the nitrogen atom obtained with GUGA-FCIQMC and CCSD(T) Knowles, Hampel, and Werner ( 1993 , 2000 ) for different basis set sizes cc-pV n 𝑛 n Z, n = 2 , 3 , 4 , 5 , 6 𝑛 2 3 4 5 6 n=2,3,4,5,6 (2 = D, 3 = T, 4 = Q) and CBS limit extrapolations, with Eq. ( 65 ) using the n = 3 , 4 , 5 𝑛 3 4 5 n=3,4,5 and 6 6 6 results, compared with experimental results Kramida  et al. ( 2018 ); Gallagher and Moore ( 1993 ); Eriksson and Pettersson ( 1971 ) . Energies are given in atomic units.

	${}_{}{}^{2}D_{}^o-{}_{}{}^{4}S_{}^o$ spin gap			N+ ${}_{}{}^{3}P_0^{}$ - N ${}_{}{}^{4}S_{}^o$ IP
$n$	CCSD(T)	GUGA-FCIQMC		CCSD(T)	GUGA-FCIQMC
2	0.1061699	0.099951(11)0		0.5216483	0.5215711(23)
3	0.1013301	0.0923719(90)		0.5310023	0.5310210(83)
4	0.0994312	0.0896661(73)		0.5333069	0.5333855(25)
5	0.0986498	0.0885878(67)		0.5341573	0.534204(12)0
6	0.0983705	0.0881762(69)		0.5344735	0.5345070(62)
CBS	()00.097950(35)	0.0875830(80)		(1)0.534987(43)	0.534971(13)0
Experiment	0.08746(37)00			0.5341192(15)
$\Delta E$	()0.01034(40)	0.00003(38)00		(1)0.000868(46)	0.000852(16)0

Table 9. Table 9: N 2 dissociation energy obtained with GUGA-FCIQMC in the frozen-core approximation in a cc-pV n 𝑛 n Z basis set and all-electron calculations in a cc-pCV n 𝑛 n Z basis set for increasing cardinal number n 𝑛 n compared to CCSD(T) Werner  et al. ( 2012 , 2015 ); Knowles, Hampel, and Werner ( 1993 , 2000 ) and experimental results Huber and Herzberg ( 1979 ) . The frozen-core experimental results are corrected for scalar relativistic, spin-orbit and frozen-core effects according to ( 87 ; 88 ) . The CBS limit is obtained with Eq. ( 65 ) for the n 𝑛 n = 3, 4 and 5 data points where available and with a Helgaker two-point extrapolation Helgaker  et al. ( 1997 ) for the T and Q GUGA-FCIQMC all-electron results. All energies are given in E h subscript 𝐸 ℎ E_{h} .

	Frozen-core cc-pV$n$Z			All-electron cc-pCV$n$Z
$n$	CCSD(T)	GUGA-FCIQMC		CCSD(T)	GUGA-FCIQMC
2	0.3184752()	0.3198257(42)		0.3203232()	0.3215439(80)
3	0.3448342()	0.345412(26)0		0.3472842()	0.347445(35)0
4	0.3551440()	0.355565(55)0		0.3566507()	0.356759(46)0
5	0.3587423()	0.358984(49)0		0.3601695()
CBS	00.362603(47)a	0.362797(53)a0		0.3634857a()	0.363555(84)b0
Experiment	0.362700(10)c			0.364002(10)d
$\Delta E$	00.000097(57)	-0.000097(63)0-		00.000516(10)	0.000447(94)0

Table 10. Table 10: Results of the least-squares polynomial fit f ​ ( n ) = a + b ⋅ n c 𝑓 𝑛 𝑎 ⋅ 𝑏 superscript 𝑛 𝑐 f(n)=a+b\cdot n^{c} of the timestep Δ ​ τ ​ ( n ) Δ 𝜏 𝑛 \Delta\tau(n) and time per iteration t i ​ t ​ e ​ r ​ ( n ) subscript 𝑡 𝑖 𝑡 𝑒 𝑟 𝑛 t_{iter}(n) as a function of the cardinal number n 𝑛 n of the basis set for the Slater determinant (SD) based and spin-adapted (CSF) FCIQMC implementation.

		$a$	$b$	$c$
$\Delta \tau$	SD:	$1.25\cdot {10}^{-4}$	$3.53\cdot {10}^{-2}$	-2.66
	CSF:	$1.05\cdot {10}^{-5}$	$3.95\cdot {10}^{-3}$	-2.00
$t_{iter}$	SD:	$1.14\cdot {10}^{-7}$	$5.04\cdot {10}^{-9}$	-2.74
	CSF:	$1.75\cdot {10}^{-7}$	$2.08\cdot {10}^{-9}$	-3.41

Table 11. Table 11: Averaged timestep Δ ​ τ Δ 𝜏 \Delta\tau and time per iteration t 𝑡 t ratios of CSF- and SD-based FCIQMC calculations for the nitrogen atom and dimer with sample sizes n s subscript 𝑛 𝑠 n_{s} . The standard errors δ ​ t 𝛿 𝑡 \delta t and δ ​ Δ ​ τ 𝛿 Δ 𝜏 \delta\Delta\tau are also given.

System	$n_s$	$\Delta {\tau }_{SD}/\Delta {\tau }_{CSF}$	$\delta \Delta \tau$	$t_{SD}/t_{CSF}$	$\delta t$
N	10	2.55	0.12	0.90	0.08
N2	12	3.87	0.48	0.78	0.05

Table 12. Table 12: GUGA-FCIQMC, CASPT2 Aquilante  et al. ( 2016 ) and CCSD(T) Werner  et al. ( 2012 , 2015 ) results for the F 2 − F 4 superscript 𝐹 2 superscript 𝐹 4 {}^{2}F-{{}^{4}F} spin gap of Co in an ANO Aquilante  et al. ( 2016 ) and cc-basis set Werner  et al. ( 2012 , 2015 ) compared with the experimental values Kramida  et al. ( 2018 ); Sugar and Corliss ( 1985 ); Pickering and Thorne ( 1996 ) . CBS limit extrapolations were obtained with Eq. ( 65 ) with the used data points of the basis sets in parentheses. The CCSD(T) results are obtained by running MOLPRO in default mode, without further specification of the reference configuration. Further investigations are being conducted on the performance of CCSD(T) by varying the reference determinant, and will be reported elsewhere.

		Co ${}_{}{}^{2}F-{}_{}{}^{4}F$ spin gap $|E_h$
	Basis set	GUGA-FCIQMC	CASPT2	CCSD(T)
ANO-basis	VDZP	0.04895(32)0	0.04667
	VTZP	0.043358(40)	0.04373
	VQZP	0.03655(29)0	0.03675
	Full	0.03626(21)0	0.03565
	Primitive		0.03565
	CBS	0.03158(50)a0	0.03165a
cc-basis	n = D	0.046448(88)		0.1057278(00)
	n = T	0.03855(22)0		0.1054354(00)
	n = Q	0.03685(27)0		0.1052032(00)
	CBS	0.0347(50)b00		0.1050555b(00)
	Experiment		0.032285
	$\Delta E_{ANO}$	0.00070(50)0	0.000630
	$\Delta E_{cc}$	-0.00242(50)0-		0.0729115(92)

Table 13. Table 13: Relevant terms of the exchange contribution to diagonal matrix elements.

$d_i|d_j$		0	1	2	3
0		0	0	0	0
1		0	$-\frac{1}{2}\left(1+A(b_i,2,0)A(b_j,-1,1){\prod }_{k}f(b_k,d_k)\right)$	$-\frac{1}{2}\left(1-A(b_i,2,0)A(b_j,3,1){\prod }_{k}f(b_k,d_k)\right)$	-1
2		0	$-\frac{1}{2}\left(1-A(b_i,0,2)A(b_j,-1,1){\prod }_{k}f(b_k,d_k)\right)$	$-\frac{1}{2}\left(1+A(b_i,0,2)A(b_j,3,1){\prod }_{k}f(b_k,d_k)\right)$	-1
3		0	-1	-1	-2

Table 14. Table 14: Modified matrix element contributions Δ ¯ ​ ( i ) ¯ Δ 𝑖 \underline{\Delta}(i) and Δ ¯ ​ ( j ) ¯ Δ 𝑗 \overline{\Delta}(j) necessary for the on-the-fly matrix element calculation during the excitation process in the FCIQMC algorithm.

$d^{\prime }$	$d$	$\underset{¯}{R}$	$R\overline{L}$	$\underset{¯}{\Delta }\left(i\right)$		$\overline{L}$	$\underset{¯}{R}L$	$\overline{\Delta }\left(j\right)$
1	0	1	-$tA(0,2)$	-$tA(0,2)$		1	$tA(2,0)$	$tA(2,0)$
2	0	1	$tA(2,0)$	$tA(2,0)$		1	-$tA(0,2)$	-$tA(0,2)$
3	1	$A(1,0)$	-$tA(-1,0)$	$tA(-1,1)$		$A(0,1)$	-$tA(2,1)$	-$tA(2,0)$
3	2	$A(1,2)$	-$tA(3,2)$	-$tA(3,1)$		$A(2,1)$	$tA(0,1)$	tA(0,2)
$d^{\prime }$	$d$	$\underset{¯}{L}$	$\overline{R}L$	$\underset{¯}{\Delta }\left(i\right)$		$\overline{R}$	$R\underset{¯}{L}$	$\overline{\Delta }\left(j\right)$
0	1	1	-t$A(-1,1)$	-$tA(-1,1)$		1	$tA(2,0)$	$tA(2,0)$
0	2	1	$tA(3,1)$	$tA(3,1)$		1	-$tA(0,2)$	-$tA(0,2)$
1	3	$A(2,1)$	$tA(0,1)$	$tA(0,2)$		$A(0,1)$	-$tA(2,1)$	-$tA(2,0)$
2	3	$A(0,1)$	-$tA(2,1)$	-$tA(2,0)$		$A(2,1)$	$tA(0,1)$	$tA(0,2)$

Table 15. Table 15: Modified two-body terms at the single overlap site k 𝑘 k used to formulate the two-body contribution to the single excitations multiplicatively and allow an on-the-fly matrix element calculation thereof.

	Usual value					Modified Value					Multiplicative Factor $r_k$
	$R$		$L$			$\underset{¯}{R}\overline{R}$		$\underset{¯}{L}\overline{L}$					$i>j$
$d^{\prime }d|\Delta b$	-1	+1	-1	+1		- 1	+1	-1	+1		-1	+1	-1	+1
00	-1	-1	-1	-1		0	0	0	0		-0	-0	-0	-0
11	-1	$C(b,0)$	$C(b,1)$	-1		1	0	0	1		-1	-0	-0	-1
12	-$1/\left(b+2\right)$	$-$	$1/\left(b+1\right)$	$-$		1	$-$	1	$-$		-$\left(b+2\right)$	$-$	$\left(b+1\right)$	$-$
21	$-$	$1/b$	$-$	-$1/\left(b+1\right)$		$-$	1	$-$	1		$-$	-$b$	$-$	-$\left(b+1\right)$
22	$C(b,2)$	-1	-1	$C(b,1)$		0	1	1	0		-0	-1	-1	-0
33	-1	-1	-1	-1		1	1	1	1		-1	-1	-1	-1

Table 16. Table 16: L ​ L ¯ , L ¯ ​ L ∗ 𝐿 ¯ 𝐿 ¯ 𝐿 superscript 𝐿 L\underline{L},\underline{L}L^{*} and R ¯ ​ R , R ​ R ¯ ∗ ¯ 𝑅 𝑅 𝑅 superscript ¯ 𝑅 \underline{R}R,R\underline{R}^{*} contributions at the overlap range start. ∗ indicates the sign change of the x = 1 𝑥 1 x=1 matrix element, depending on the order of operators. N k ′ = N k − 2 superscript subscript 𝑁 𝑘 ′ subscript 𝑁 𝑘 2 N_{k}^{\prime}=N_{k}-2 for two lowering generators and N k ′ = N k + 2 subscript superscript 𝑁 ′ 𝑘 subscript 𝑁 𝑘 2 N^{\prime}_{k}=N_{k}+2 for two raising generators. And R ​ R 𝑅 𝑅 RR and L ​ L 𝐿 𝐿 LL intermediate segments in the overlap region of a double excitation, depending on Δ ​ b k − 1 Δ subscript 𝑏 𝑘 1 \Delta b_{k-1} .

$L\underset{¯}{L}/\underset{¯}{L}{L}^{\ast }$					$\underset{¯}{R}R/R{\underset{¯}{R}}^{\ast }$					$RR/LL$
$\Delta b_{k-1}:$		-1	+1				-1	+1				0	-2	+2
$d^{\prime }$	$d$	$\Delta b_k$			$d^{\prime }$	$d$	$\Delta b_k$			$d^{\prime }$	$d$	$\Delta b_k$
0	1	0	+2		1	0	-2	0		0	0	0	-2	+2
0	2	-2	0		2	0	0a	+2b		1	1	0	-2	+2
1	3	-2	0		3	1	0	+2		2	1	+2b	0a	$-$
2	3	0a	+2b		3	2	-2	0		1	2	-2	$-$	+2b
										2	2	0	-2a	+2b
										3	3	0	-2	+2

Table 17. Table 17: Segment value restriction for the end of the overlap range for two lowering generators L ¯ ​ L , L ​ L ¯ ∗ ¯ 𝐿 𝐿 𝐿 superscript ¯ 𝐿 \overline{L}L,L\overline{L}^{*} , and two raising generators R ​ R ¯ , R ¯ ​ R ∗ 𝑅 ¯ 𝑅 ¯ 𝑅 superscript 𝑅 R\overline{R},\overline{R}R^{*} , depending on Δ ​ b k − 1 Δ subscript 𝑏 𝑘 1 \Delta b_{k-1} value. N k ′ = N k ± 1 superscript subscript 𝑁 𝑘 ′ plus-or-minus subscript 𝑁 𝑘 1 N_{k}^{\prime}=N_{k}\pm 1 depending on the generator type. ∗ indicated that the x = 1 𝑥 1 x=1 matrix element contribution has opposite sign for exchanged order of generators.

$\overline{L}L/L{\overline{L}}^{\ast }$					$R\overline{R}/\overline{R}{R}^{\ast }$
$d^{\prime }$	$d$	$\Delta b_{k-1}$	$\Delta b_k$		$d^{\prime }$	$d$	$\Delta b_{k-1}$	$\Delta b_k$
1	0	0	-1		0	1	0	+1
		+2	+1				-2	-1
2	0	0	+1a		0	2	0	-1
		-2	-1				+2	+1
3	1	0	+1		1	3	0	-1
		-2	-1				+2	+1
3	2	0	-1		2	3	0	+1a
		+2	+1				-2	-1

Table 18. Table 18: R ​ L ¯ ¯ 𝑅 𝐿 \underline{RL} starting segments, where there is no change in the orbital occupation number N k ′ = N k superscript subscript 𝑁 𝑘 ′ subscript 𝑁 𝑘 N_{k}^{\prime}=N_{k} . And R ​ L ¯ ¯ 𝑅 𝐿 \overline{RL} end segment restrictions, depending on the Δ ​ b k − 1 Δ subscript 𝑏 𝑘 1 \Delta b_{k-1} value.

		$\underset{¯}{R}\underset{¯}{L}$		$\overline{R}\overline{L}$
$d^{\prime }$	$d$	$\Delta b_k$		$\Delta b_{k-1}$
1	1	0		0
2	1	+2a		-2
1	2	-2		+2
2	2	0		0
3	3	0b		0b

Equations

$$[\hat{H},{\hat{\mathbf{S}}}^2]=0,$$

$$\hat{H}=\sum _{ij}^n{t}_{ij}\sum _{\sigma =\uparrow ,\downarrow }{a}_{i\sigma }^{†}{a}_{j\sigma }+\frac{1}{2}\sum _{ijkl}^n⟨ik|r_{12}^{-1}|jl⟩\sum _{\sigma ,\tau =\uparrow ,\downarrow }{a}_{i\sigma }^{†}{a}_{k\tau }^{†}{a}_{l\tau }{a}_{j\sigma },$$

$$a_{i\sigma }^{†}{a}_{k\tau }^{†}{a}_{l\tau }{a}_{j\sigma }=a_{i\sigma }^{†}{a}_{j\sigma }{a}_{k\tau }^{†}{a}_{l\tau }-{\delta }_{jk}{\delta }_{\sigma \tau }{a}_{i\sigma }^{†}{a}_{l\sigma },$$

$$\sum _{\sigma }{a}_{i\sigma }^{†}{a}_{j\sigma }={\hat{E}}_{ij}$$

$$\sum _{\sigma \tau }{a}_{i\sigma }^{†}{a}_{k\tau }^{†}{a}_{l\tau }{a}_{j\sigma }={\hat{E}}_{ij}{\hat{E}}_{kl}-{\delta }_{jk}{\hat{E}}_{il}={\hat{e}}_{ij,kl}.$$

$$\hat{H}=\sum _{ij}{t}_{ij}{\hat{E}}_{ij}+\frac{1}{2}\sum _{ij,kl}{V}_{ij,kl}{\hat{e}}_{ij,kl}.$$

$$[{\hat{E}}_{ij},{\hat{E}}_{kl}]=$$

$$=$$

$$[{\hat{E}}_{ij},{\hat{E}}_{kl}]=$$

$$U(n)\supset U(n-1)\supset \cdots \supset U(2)\supset U(1),$$

$$[\mathbf{m}]=\left[\begin{array}{ccccccc}{m}_{1,n}& & m_{2,n}& \cdots & m_{n-1,n}& & m_{n,n}\\ & m_{1,n-1}& & \cdots & & m_{n-1,n-1}& \\ & \ddots & & \cdots & & \text{iddots}& \\ & & m_{1,2}& & m_{2,2}& & \\ & & & m_{1,1}& & & \end{array}\right].$$

$$m_{i,j+1}\ge m_{ij}\ge m_{i+1,j+1},$$

$${\hat{A}}_{i\sigma ,j\tau }={\hat{a}}_{i\sigma }^{†}{\hat{a}}_{j\tau };i,j=1,\dots ,n;\sigma ,\tau =\uparrow ,\downarrow$$

$$[{\hat{A}}_{i\sigma ,j\tau },{\hat{A}}_{i^{\prime }{\sigma }^{\prime },j^{\prime }{\tau }^{\prime }}]={\delta }_{j{i}^{\prime }}{\delta }_{\tau {\sigma }^{\prime }}{\hat{A}}_{i\sigma ,j^{\prime }{\tau }^{\prime }}-{\delta }_{i{j}^{\prime }}{\delta }_{\sigma {\tau }^{\prime }}{\hat{A}}_{i^{\prime }{\sigma }^{\prime },j\tau }.$$

$${\hat{E}}_{ij}=\sum _{\sigma =\uparrow ,\downarrow }{\hat{A}}_{i\sigma ,j\sigma }\text{and}{\hat{\mathcal{E}}}_{\sigma \tau }=\sum _{i=1}^n{\hat{A}}_{i\sigma ,i\tau }$$

$$U(2n)\supset U(n)\otimes U(2),$$

$$\sum _{i=1}^n{m}_{in}=N\text{and}\frac{1}{2}\sum _{i=1}^n{\delta }_{1,m_{in}}=S,$$

$${\footnotesize\begin{bmatrix}\phantom{0}\quad\phantom{0}\quad\phantom{0}\quad\phantom{0}\quad\phantom{0}\quad\phantom{0}\quad\phantom{0}\quad\phantom{0}\\ 2\quad 2\quad 1\quad 1\quad 0\quad 0\quad 0\quad 0\\ 2\quad 2\quad 1\quad 1\quad 0\quad 0\quad 0\\ 2\quad 1\quad 1\quad 0\quad 0\quad 0\\ 2\quad 1\quad 1\quad 0\quad 0\\ 2\quad 1\quad 0\quad 0\\ 1\quad 1\quad 0\\ 1\quad 0\\ 1\end{bmatrix}\equiv\left[\begin{array}[]{ccc}a_{i}&b_{i}&c_{i}\\ \hline\cr 2&2&4\\ 2&2&3\\ 1&2&3\\ 1&2&2\\ 1&1&2\\ 0&2&1\\ 0&1&1\\ 0&1&0\end{array}\right]\rightarrow\left[\begin{array}[]{ccc}\Delta a_{i}&\Delta b_{i}&\Delta c_{i}\\ \hline\cr 0&\phantom{-}0&1\\ 1&\phantom{-}0&0\\ 0&\phantom{-}0&1\\ 0&\phantom{-}1&0\\ 1&-1&1\\ 0&\phantom{-}1&0\\ 0&\phantom{-}0&1\\ 0&\phantom{-}1&0\end{array}\right]}$$

$$a_i+b_i+c_i=i,(i=1,\dots ,n)$$

$$a=a_n=\frac{1}{2}N-S,b=b_n=2S,c=c_n=n-a-b=n-\frac{1}{2}N-S,$$

$$N_{CSF}=\frac{b+1}{n+1}\left(\genfrac{}{}{0ex}{}{n+1}{a}\right)\left(\genfrac{}{}{0ex}{}{n+1}{c}\right)=\frac{2S+1}{n+1}\left(\genfrac{}{}{0ex}{}{n+1}{\frac{N}{2}-S}\right)\left(\genfrac{}{}{0ex}{}{n+1}{n-\frac{N}{2}-S}\right).$$

$$\frac{N_{SD}}{N_{CSF}}\approx \frac{\sqrt{\pi n}n}{2(2S+1)},$$

$$N_i=2a_i+b_i,S_i=\frac{1}{2}{b}_i.$$

$$N_{DRT}=$$

$$=$$

$$⟨m^{\prime }|\hat{H}|m⟩=\sum _{ij}{t}_{ij}⟨m^{\prime }|{\hat{E}}_{ij}|m⟩+\frac{1}{2}\sum _{ij,kl}{V}_{ij,kl}⟨m^{\prime }|{\hat{e}}_{ij,kl}|m⟩.$$

$${\hat{E}}_{ij}|m⟩=\sum _{m^{\prime }}|m^{\prime }⟩⟨m^{\prime }|{\hat{E}}_{ij}|m⟩,$$

$$N_k^{\prime }=N_k\pm 1\text{and}{b}_k^{\prime }=b_k\pm 1\iff S_k^{\prime }=S_k\pm \frac{1}{2}\text{for}k\in S_0,$$

$$a_k^{\prime }=a_k,b_k^{\prime }=b_k+1,c_k^{\prime }=c_k-1,$$

$$a_k^{\prime }=a_k+1,b_k^{\prime }=b_k-1,c_k^{\prime }=c_k$$

$$a_k^{\prime }=a_k-1,b_k^{\prime }=b_k+1,c_k^{\prime }=c_k$$

$$a_k^{\prime }=a_k,b_k^{\prime }=b_k-1,c_k^{\prime }=c_k+1$$

Full text

Quantum Monte Carlo, SU(2) symmetry

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, $L_{z}$ angular momentum and $S_{z}$ projected spin conservation.

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

[TABLE]

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

In addition to the above-mentioned Hilbert space size reduction and conservation of the total spin $S$, solving for the eigenstates of $\hat{H}$ in a simultaneous spin-eigenbasis of $\hat{\mathbf{S}}^{2}$ 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 $\hat{\mathbf{S}}^{2}$ 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 $S$ 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 and excited states wavefunctions are formulated nonlinearly based on the DRT, conserving the total spin $S$.

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 N2 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 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 $S$ 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 $\hat{H}$ 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

[TABLE]

where $t_{ij}=\langle i|\hat{h}|j\rangle$.With the reformulation

[TABLE]

we can define

[TABLE]

and

[TABLE]

as the singlet one- and two-body excitation operators Helgaker, Jørgensen, and Olsen (2000), which do not change $S$ and $m_{s}$ upon acting on a state, $\ket{S,m_{s}}$, 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)

[TABLE]

where $V_{ij,kl}=\langle ik|r_{12}^{-1}|jl\rangle$. 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 $U(n)$ Paldus (1974, 1975, 1976), $n$ being the number of spatial orbitals. The commutator of the spin-preserving excitation operators $\hat{E}_{ij}$ can be calculated as

[TABLE]

which is the same as for the basic matrix units and the generators of the unitary group $U(n)$.

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 $\hat{H}$ does not commute with the generators of $U(n)$, 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).

We only want to recap the most important concepts of the UGA here and refer the interested reader to the pioneering work of Paldus (53) and Shavitt (42; 43; 69).

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 $U(n)$. The group $U(n)$ has $n^{2}$ generators, $E_{ij}$, and a total of $n$ Casimir operators, commuting with all generators of the group, and the GT basis is based on the group chain

[TABLE]

where $U(1)$ is Abelian and has one-dimensional irreducible representations (irreps). Each subgroup

$U(n-1),U(n-2),\dots,U(1)$ has $n-1,n-2,\dots,1$ Casimir operators, resulting in a total of $n(n+1)/2$ 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 $n(n+1)/2$ integers related to the eigenvalues of the invariants. Thus, based on the branching law of Weyl Weyl (1931, 1946), a general $N$-electron CSF can be represented by a Gel’fand pattern Gel’fand and Cetlin (1950a)

[TABLE]

The integers $m_{ij}$ in the top row (and all subsequent rows) of (8) are nonincreasing, $m_{1n}\geq m_{2n}\geq\dots\geq m_{nn}$, and the integers in the subsequent rows fulfill the condition

[TABLE]

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

The $n$ non-increasing integers of the top row of Eq. (8), $\mathbf{m}_{n}=\left(m_{1n},m_{2n},\dots,m_{nn}\right)$, are called the highest weight or *weight vector * of the representation and specify the chosen irrep of $U(n)$; the following $n-1$ rows uniquely label the states belonging to the chosen irrep.

In CI calculations one usually employs a one-particle basis of $2n$ spin-orbitals with creation $\hat{a}_{i\sigma}^{\dagger}$ and annihilation $\hat{a}_{j\tau}$ operators of electrons in spatial orbital $i,j$ with spin $\sigma,\tau$. The $(2n)^{2}$ operators

[TABLE]

can be associated with the generators of $U(2n)$ with the commutation relation

[TABLE]

The partial sums over spin or orbital indices of these operators

[TABLE]

are related to the orbital $U(n)$ and spin $U(2)$ generators. Since we deal with fermions we have to restrict ourselves to the totally antisymmetric representations of $U(2n)$, denoted as $\Gamma\{1^{2n}\}$. Since the molecular Hamiltonian (5) is spin independent, we can consider the proper subgroup of the direct product of the spin-free orbital space $U(n)$, with $n^{2}$ generators $E_{ij}$, and the pure spin space $U(2)$ with the four generators ${\mathcal{E}}_{\sigma\tau}$ Paldus (1974), given as

[TABLE]

where the representations of $U(n)$ and $U(2)$ are mutually conjugate Paldus (2006); Matsen (1974, 1964); Moshinsky (1968).

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

II.2 The Paldus tableau

The consequence of the mutually conjugate relationship between $U(n)$ and $U(2)$ irreps for electronic structure calculations is that the integers $m_{ij}$ in a Gel’fand pattern (8) for $U(n)$ are related to occupation numbers of spatial orbitals. This means they are restricted to $0\leq m_{ij}\leq 2$, due to the Pauli exclusion principle. The highest weight, $\mathbf{m}_{n}$, indicates the chosen electronic state with the conditions

[TABLE]

with $N$ being the total number of electrons and the number of singly occupied orbitals, $\delta_{1,m_{ij}}$ is equal to twice the total spin value $S$.

This insight led Paldus Paldus (1974) to the more compact formulation of a GT state by a table of $3n$ integers. It is sufficient to count the appearances $2^{\prime}s$, $1^{\prime}s$ and $0^{\prime}s$ in each row $i$ of a Gel’fand pattern and store this information, denoted by $a_{i},b_{i}$ and $c_{i}$ in a table, named a Paldus tableau.

The first column, $a_{i}$, contains the number of doubly occupied orbitals, the second column, $b_{i}$, the number of singly occupied and the last one, $c_{i}$, the number of empty orbitals, as shown by the example of an $n=8,N=6,S=1$ state:

[TABLE]

where the differences $\Delta x_{i}=x_{i}-x_{i-1}$, with $x=a,b,c$, of subsequent rows are also indicated. For each row the condition

[TABLE]

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

[TABLE]

completely specifying the chosen electronic state as an irrep of $U(n)$.

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

[TABLE]

As it can be seen from Eq. (18), 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 $N=n$ can be estimated by Stirling’s formula (for sufficiently large $n$ and $N$) as

[TABLE]

which shows orbital dependent, $\sim n^{3/2}$, decrease of the efficient Hilbert space size for a spin-adapted basis. The Paldus tableau also emphasizes the cumulative aspects of the coupling between electrons, with the i-th row providing information on number of electrons, $N_{i}$ (up to i-th level) and the spin, $S_{i}$, by

[TABLE]

As can be seen in Eq. 15, there are four permissible difference vectors $[\Delta a_{i},\Delta b_{i},\Delta c_{i}]$ ($\Delta x_{i}=x_{i}-x_{i-1}$, with $x=a,b,c$) between consecutive rows of a Paldus tableau, 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” [math]’th row $i=0$, an empty spatial orbital is indicated by $d_{i}=0$, a “positively spin-coupled” orbital, $\Delta S_{i}=1/2$, by $d_{i}=1$, a “negatively spin-coupled”, $\Delta S_{i}=-1/2$, by $d_{i}=2$ and a doubly occupied spatial orbital by $d_{i}=3$. To retain physically allowed states the condition $S_{i}\geq 0,\,\forall i$ applies. (As a side note: Another common notation —e.g. in Molcas—is to indicate positive spin-coupling as $d_{i}=u$, negative spin-coupling by $d_{i}=d$ and a doubly occupied orbital by $d_{i}=2$.)

The step-value $d_{i}$ in Tab. 1 is given by $d_{i}=3\Delta a_{i}+\Delta b_{i}$ and the collection of all $d_{i}$ into the step-vector $\mathbf{d}$ 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 $U(n)$ by constructing all possible distinct step-vectors $\ket{\mathbf{d}}$ which lead to the same top-row of the Paldus tableau (17), specifying the chosen irrep with definite spin and number of electrons, with the restriction $S_{i}\geq 0,\,\forall i$.

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 tableaux specifying the CSFs of a chosen irrep of $U(n)$. Instead of all possible Paldus tableaux, 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)

[TABLE]

with $d=\min(a,c)=\min(N/2-S,n-N/2-S)$, which is drastically smaller than the total number of possible CSFs (18) or Slater determinants (without any symmetry restrictions) as seen in Fig. 1. Each row in the DRT is identified by a pair of indices $(i,j)$, with $i=a_{ij}+b_{ij}+c_{ij}$ being the level index, related to the orbital index and $j$ being the lexical row index such that $j<j^{\prime}$ if $a_{ij}>a^{\prime}_{ij}$ or if $a_{ij}=a^{\prime}_{ij}$ and $b_{ij}>b^{\prime}_{ij}$.

A simple example of the DRT of a system with $n=3$, $N=4$ and $S=0$ is shown in Table 2.

Relations between elements of the DRT belonging to two neighboring levels $k$ and $k-1$ are indicated by the so called downward, $k_{d_{k}}$, and upward, $l_{d_{k}}$, chaining indices, with $d_{k}=0,1,2,3$. These indices indicate the connection to a lexical row index in a neighboring level by a step-value $d_{k}$, 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.

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 $j$, starting at the unique head node at the top, which corresponds to the highest row $(a,b,c)$. It ends at the second unique null row $(0,0,0)$, which is called the tail of the graph. Vertices with the same $i$-value of Table 2 are at the same level on this grid. The highest $i$-value is on top and the lowest at the bottom. Vertices also have left-right order with respect to their $a_{i}$ value and vertices that share the same $a_{i}$ value are further ordered—still horizontally—with respect to their $b_{i}$ value. With the above mentioned ordering of the vertices according to their $a_{i}$ and $b_{i}$ values, the slope of each arc is in direct correspondence to the step-value $d_{i}$, connecting two vertices. $d=0$ corresponds to vertical lines, and the tilt of the other arcs increases with the step-value $d_{i}$.

Each CSFs in the chosen irrep of $U(n)$, 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 $\ket{\mathbf{d}_{1}}=\ket{3,3,0}$ and $\ket{\mathbf{d}_{2}}=\ket{3,1,2}$ in step-vector representation. Such a walk spans $n$ arcs (number of orbitals) and visits one node at each level $i$. There is a direct correspondence between the Paldus tableau, Gel’fand patterns and directed walks on Shavitt graphs for representing all possible CSFs in a chosen irrep of $U(n)$.

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, $\ket{m^{\prime}}$ and $\ket{m}$, is given by:

[TABLE]

The matrix elements, $\langle{m^{\prime}}|{\hat{E}_{ij}}|{m}\rangle$ and $\langle{m^{\prime}}|{\hat{e}_{ij,kl}}|{m}\rangle$, provide the coupling coefficients between two given CSFs and $t_{ij}$ and $V_{ij,kl}$ are the integral contributions. The coupling coefficients are independent of the orbital shape and only depend on the involved CSFs, $\ket{m^{\prime}}$ and $\ket{m}$- 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) $\hat{E}_{ij}$, between two GT states (CSFs), $\langle{m^{\prime}}|{\hat{E}_{ij}}|{m}\rangle$. The generators are classified according to their indices, with $\hat{E}_{ii}$ being diagonal weight (W) and $\hat{E}_{ij}$ with $i<j$ being raising (R) and $i>j$ lowering (L) operators (or generators). In contrast to Slater determinants, $\hat{E}_{ij}$ applied to $\ket{m}$ yields a linear combination of CSFs $\ket{m^{\prime}}$,

[TABLE]

with an electron moved from spatial orbital $j$ to orbital $i$ without changing the spin of the resulting states $\ket{m^{\prime}}$. They are called raising (lowering) operators since the resulting $\ket{m^{\prime}}$ will have a higher (lower) lexical order than the starting CSF $\ket{m}$.

The distance, $S_{0}$, from $\min(i,j)-1$ to $\max(i,j)$, is an important quantity and is called the range of the generator $\hat{E}_{ij}$. For the one-body term in (5) Shavitt Shavitt (1977) was able to show that the walks on the graph, representing the CSFs $\ket{m}$ and $\ket{m^{\prime}}$, must coincide outside of this range $S_{0}$ to yield a non-zero matrix element. The two vertices in the DRT graph, related to orbital $i-1$ and $j$ (with $i<j$) represent the points of separation of the walks and they are named loop head and loop tail. And the matrix element $\langle{m^{\prime}}|{\hat{E}_{ij}}|{m}\rangle$ only depends on the shape of the loop formed by the two graphs in the range $S_{0}$, shown in Fig. 3.

Shavitt Shavitt (1978) showed that the relations

[TABLE]

between $\ket{m}$ and $\ket{m^{\prime}}$ must be fulfilled to yield a nonzero matrix element ($N^{\prime}_{k}=N_{k}+1$ for a raising and $N^{\prime}_{k}=N_{k}-1$ 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:

[TABLE]

where $\Delta b_{k}=b_{k}-b_{k}^{\prime}$ and for lowering generators L:

[TABLE]

At each vertex of the loop in range $k$ one of the relations (25-28) 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 $\hat{E}_{ij}$ can be factorized in a product, where each term corresponds to a segment of the loop in the range $S_{0}$ and is given by

[TABLE]

where $b_{k}$ is the $b$ value of state $\ket{m}$ at level $k$. $W(Q_{k};d^{\prime}_{k},d_{k},\Delta b_{k},b_{k})$ additionally depends on the segment shape of the loop at level $k$, determined by the type of the generator $Q_{k}=W,R,L$, the step values $d^{\prime}_{k}$ and $d_{k}$ and $\Delta b_{k}=b_{k}-b^{\prime}_{k}$. 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 $\hat{E}_{ij}$—an over/under-bar indicates the loop head/tail—depending on the segment shape symbol, the step-values and the $b$-value are given in terms of the auxiliary functions

[TABLE]

III.2 Two-Body Matrix Elements

The matrix elements of the two-body operators $\hat{e}_{ij,kl}$ are more involved than the one-body operators, especially the product of singlet excitation generators, $\hat{E}_{ij}\hat{E}_{kl}$. Similar to the one-electron operators, the GT states $\ket{m}$ and $\ket{m^{\prime}}$ must coincide outside the total range $\min(i,j,k,l)$ to $\max(i,j,k,l)$ for $\langle{m^{\prime}}|{\hat{e}_{ij,kl}}|{m}\rangle$ to be nonzero. The form of the matrix element depends on the overlap range of the two ranges

[TABLE]

One possibility to calculate the matrix element would be to sum over all possible intermediate states, $\ket{m^{\prime\prime}}$,

[TABLE]

but in practice this is very inefficient. For non-overlapping ranges $S_{1}=\emptyset$ the matrix element just reduces to the product

[TABLE]

where $\ket{m^{\prime\prime}}$ must coincide with $\ket{m}$ in the range $(i,j)$ and with $\ket{m^{\prime}}$ in range $(k,l)$. 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 $S_{1}\neq\emptyset$, we define the non-overlap range

[TABLE]

where the same restrictions and matrix elements as for one-body operators apply. In the overlap range, $S_{1}$, different restrictions for the visited Paldus tableau vertices $p$ 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

[TABLE]

For two lowering generators (LL):

[TABLE]

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

[TABLE]

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 (29)

[TABLE]

where $S_{1}$ and $S_{2}$ are the overlap (31) and non-overlap (34) ranges defined above.

$W(Q_{p};d^{\prime}_{p},d_{p},\Delta b_{p},b_{p})$ are the already defined single operator segment values, listed in Table 3, and $W_{x}(Q_{p};d^{\prime}_{p},d_{p},\Delta b_{p},b_{p})$ are new segment values of the overlap range (their listing is omitted for brevity here, but can be found in Refs. [(69; 74)]. The sum over two products in $S_{1}$ corresponds to the singlet coupled intermediate states ($x=0$), with a nonzero contribution if $\Delta b_{p}=0,\forall p\in S_{1}$ and the triplet intermediate coupling ($x=1$).

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. [(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. [(1; 2)] and [(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. (i)

Efficient storage of the spin-adapted basis 2. (ii)

Efficient excitation identification and matrix element computation 3. (iii)

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. [(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 $\hat{H}$ than for SDs, it is possible to do so efficiently. Matrix element computation is based on the product structure of the one- (29) and two-body (44) matrix elements derived by Shavitt Shavitt (1978) explained above and presented in more detail in App. B and in Ref. [(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 $SU(2)$ 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 $S$ 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 $i$ and $j$, with probability $p(i)\,p(j|i)$. 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 $(i,\sigma)$ and an empty spin-orbital $(j,\sigma)$ is sufficient to uniquely specify the excitation $\ket{D_{j}}=a_{j,\sigma}^{\dagger}a_{i,\sigma}\ket{D_{i}}$, and to calculate the involved matrix element $\langle{D_{j}}|{\hat{H}}|{D_{i}}\rangle$. However, in a CSF basis, the choice of an occupied spatial orbital $j$, and empty or singly occupied spatial orbital $i$, only determines the type of excitation generator $\hat{E}_{ij}$ acting on an CSF basis state $\ket{m}$ as well as the involved integral contributions $t_{ij},V_{ikjk}$ and $V_{ikkj}$ of the matrix element $\langle{m^{\prime}}|{\hat{H}}|{m}\rangle$. To ensure $p(i)p(j|i)\propto\lvert t_{ij}+\frac{1}{2}\sum_{k\in\text{occ}}(V_{ikjk}-V_{ikkj})\rvert$, the occupied orbital $j$ and (partially) empty $i$ are picked in the same way as for SDs, but with an additional restriction to ensure $\hat{E}_{ij}\ket{m}\neq 0$. However, the choice of $(i,j)$ does not uniquely determine the excited CSF as there are multiple possible ones, as explained above.

As a consequence, the choice of spatial orbitals $i$ and $j$ does not determine the coupling coefficient $\langle{m^{\prime}}|{\hat{E}_{ij}}|{m}\rangle$ of the matrix element $H_{m^{\prime}m}$. Optimally, for a given $\ket{m}$ and generator $\hat{E}_{ij}$, the connected CSF $\ket{m^{\prime}}$ has to be created with a probability $p(m^{\prime}|m)$ proportional to the coupling coefficient $\langle{m^{\prime}}|{\hat{E}_{ij}}|{m}\rangle$. By ensuring $p(i)p(j|i)$ is proportional to the integral contributions and $p(m^{\prime}|m)$ to the coupling coefficients, the total spawning probability

[TABLE]

will be proportional to the magnitude of Hamiltonian matrix element $\lvert H_{m^{\prime}m}\rvert$. The efficiency of the FCIQMC algorithm depends on the ratio of the Hamiltonian matrix element $\lvert H_{m^{\prime}m}\rvert$ between two connected states and the probability $p_{s}(m^{\prime}|m)$ to choose the excitation $\ket{m}\rightarrow\ket{m^{\prime}}$, as the imaginary timestep $\Delta\tau$ of the simulation is adapted to faithfully account for all excitations

[TABLE]

In a primitive implementation, $\Delta\tau$ is determined by the “worst-case” $\max{\lvert H_{m^{\prime}m}/p_{s}(m^{\prime}|m)\rvert}$ ratio during a simulation. A less strict approach to this problem is discussed below. By choosing nonzero $\hat{E}_{ij}\ket{m}\neq 0$ and ensuring $p(m^{\prime}|m)$ 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 $\ket{m}$.

IV.2 The Branching Tree

In the spin-adapted excitation generation, after a certain generator $\hat{E}_{ij}$ is picked with a probability $p(i)p(j|i)$ based on the integral contributions of the Hamiltonian matrix element, the type of generator is determined, raising (R) if $i<j$ and lowering (L) if $i>j$. One connecting single excitation is then chosen by looping from starting orbital $\min(i,j)$ to $\max(i,j)$ and stochastically choosing a valid nonzero Shavitt graph, based on the restrictions (25-28), 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 $\underline{R}$ matrix elements. These starting segments are associated with a relative difference of the total spin $\Delta S_{i}=S_{i}(m^{\prime})-S_{i}(m)$ and $\Delta b_{i}=b_{i}(m^{\prime})-b_{i}(m)$ between the two CSFs $\ket{m}$ and $\ket{m^{\prime}}$ at level $i$, as shown in Table 4. For certain step-values ($d_{i}=0$ for raising and $d_{i}=3$ for lowering generators) two possible excited CSFs with different $\Delta b_{i}$ are possible. This can be represented pictorially as elements of a branching tree, as seen in Fig. 6 for raising generator, where the number in the boxes represent the step-value $d_{i}$ of $\ket{m}$ and the direction of the outgoing lines the $\Delta b_{i}$ value (left going lines correspond to $\Delta b_{i}=-1$ and right going ones $\Delta b_{i}=+1$). The number above the small dots represent the associated $d^{\prime}_{i}$ value of the excited $\ket{m^{\prime}}$.

The intermediate contributions to the coupling coefficients $R/L$, see Table 3, have similar properties. Depending on the current $\Delta b_{k}$ value of the excitation $\ket{m^{\prime}}$ relative to $\ket{m}$ there are branching possibilities for singly occupied spatial orbitals in $\ket{m}$, corresponding to possible spin-recouplings in the excitation range of $\hat{E}_{ij}$. An excitation with $\Delta b_{k-1}=-1$ can branch at $d_{k}=1$ values, into $d^{\prime}_{k}=1$ with $\Delta b_{k-1}=\Delta b_{k}=-1$ or change the spin-coupling to $d^{\prime}_{k}=2$ accompanied by a change to $\Delta b_{k}=+1$. At empty or doubly occupied orbitals only $d^{\prime}_{k}=d_{k}$ and $\Delta b_{k}=\Delta b_{k-1}$ leads to nonzero excitations. These relations are tabulated in Table 5 and pictorially represented in Fig. 6.

The possible single excitations of a given CSF can be represented by a branching diagram, where each node is a successive element of $d_{k}$ and a left going branch represents a $\Delta b_{k}=-1$ value and a right going branch $\Delta b_{k}=+1$. The end value $d_{j}=1$ requires an incoming $\Delta b_{j-1}=-1$ value, whereas $d_{j}=2$ requires $\Delta b_{j-1}=+1$ to ensure $\Delta b_{j}=0$ at the end of the excitation, indicated by the directions of the ingoing lines of the elements at the bottom of Fig. 6. For a raising generator both $\Delta b_{j-1}$ values are possible for $d_{j}=3$. The restrictions on the end segments $\overline{R}/\overline{L}$ are listed in Table 4 and pictorially represented in Fig. 6. These restrictions are a direct consequence of the conservation of the total spin quantum number $S$ in the GUGA.

A very simple implementation to create a single excitation $\ket{m^{\prime}}$ would be to loop from orbital $i$ to $j$ and depending on the step-value $d_{k}$ of $\ket{m}$, at each orbital $k\in(i,j)$ choose one possible $\Delta b_{k}$ path at random if there are multiple possible ones. However, this would totally neglect that there are certain branching choices which would lead to a dead end, due to incompatible end-segments $\overline{R},\overline{L}$ and would not relate the probability to create a certain CSF $\ket{m^{\prime}}$ to the magnitude of coupling coefficient.

An example of the excitation generation based on the branching tree is given in Fig. 7, for the raising generator $\hat{E}_{26}$ acting on the CSF $\ket{m}=\ket{1,0,1,2,0,1,0}$, moving an electron from spatial orbital $6$ to $2$. The left panel of Fig. 7 shows this excitation in the Shavitt graph form based on the DRT and the right panel shows the branching tree representation (with the orbitals ordered from top to bottom now, as this is the usual representation of trees.), with $\pm 1$ indicating the $\Delta b_{k}$ value associated to the possible branches. The orange path (color online) in both the Shavitt graph and branching tree representation show one valid single excitation $\ket{m^{\prime}}$ of $\hat{E}_{26}\ket{m}$. The above mentioned dead ends are indicated with dashed lines and crossed out vertices in the right panel of Fig. 7.

As one can see the number of connected CSFs $n_{m^{\prime}}$ to $\ket{m}$ via a single application of $\hat{E}_{ij}$ depends on the number of singly occupied orbitals $n_{s}$ within the excitation range $(i,j)$ and grows approximately as $n_{m^{\prime}}\approx 1.6^{n_{s}+2}$. The highest number of possible connected CSFs is given for a starting segment $\underline{R}/\underline{L}$ with two possible branches, exclusively alternating singly occupied orbitals $d_{k}=\{1,2\}$ in the excitation range with $b_{k}>0$ and an end-segment $\overline{R}/\overline{L}$ with nonzero contributions for both $\Delta b_{k}=\pm 1$. In this case the number of connected CSFs is related to the Fibonacci series and given by the Fibonacci number

[TABLE]

Calculating all possible excitations would lead to an exponential wall for highly open-shell CSFs $\ket{m}$, but since we only need to obtain one connected CSF in the excitation generation of FCIQMC, this exponential scaling is not an immediate problem. However, since the overall generation probability $\sum_{m^{\prime}}p(m^{\prime}|m)$ is normalized to unity, a specific $p(m^{\prime}|m)$ will be negligibly small for numerous possible excitations. As the timestep is directly related to this probability (46), a small $p(m^{\prime}|m)$ directly causes a lowering in the usable $\Delta\tau$ in a FCIQMC calculation.

Since this is a consequence of the inherent high connectivity of a spin-adapted basis, systems with many open-shell orbitals are difficult to treat in such a basis. In general this restricts common implementations of spin-eigenfunctions to a maximum of 18 open-shell orbitals. However, similar to the avoidance of the exponential wall associated with the FCI solution to a system, the stochastic implementation of the CSF excitation generation in FCIQMC, avoids the exponential bottle-neck caused by the high connectivity of a CSF basis.

IV.3 Remaining Switches

To avoid ending up in incompatible dead-end excitations it is convenient, for a given excitation range $(i,j)$, to determine the vector of remaining switch possibilities $s_{k}({\Delta b_{k}})$ for the $\Delta b_{k}=\pm 1$ branches. $s_{k}(\Delta b_{k})$ is the number of $d_{k^{\prime}}=1$ for $\Delta b_{k}=-1$ and $d_{k^{\prime}}=2$ for $\Delta b_{k}=+1$ to come in $k^{\prime}=k+1,\dots j-1$ (with the already mentioned restriction of $b_{k^{\prime}}>1$ for $d_{k^{\prime}}=2$ to be a valid $\Delta b_{k}=+1$ switch)

[TABLE]

The quantity $s_{k}(\Delta b_{k})$ can be used to decide if a possible $\Delta b_{k}$ branch is taken or not, depending on if it will end up in a dead-end of the branching tree.

IV.4 On-The-Fly Matrix Element Calculation

To pick the connecting CSF $\ket{m^{\prime}}$ with a probability $p(m^{\prime}|m)$ relative to the magnitude of the generator matrix element $\lvert\bra{m^{\prime}}\hat{E}_{ij}\ket{m}\rvert$ we have to investigate the matrix element $\bra{m^{\prime}}\hat{E}_{ij}\ket{m}$ between a given CSF $\ket{m}$ and an excitation $\ket{m^{\prime}}$. As the coupling coefficient is calculable as a product of terms, which depend on the type of excitation (lowering, raising) and is determined by the step-vector values $d_{k}$, $d^{\prime}_{k}$, the $b_{k}$ and the $\Delta b_{k}$ associated to each level of the excitation, see Eq. (29). One of the major advantages of the GUGA in FCIQMC is that this matrix element can be calculated on-the-fly during the creation of the excitation. As one can see in Table 3 there is a relation between the matrix element amplitude and the number of direction switches of $\Delta b_{k}$ in the excitation range. Most product contributions are of order $\mathcal{O}(1)$, except the elements related to a switch of $\Delta b_{k+1}\leftarrow-\Delta b_{k}$, which are of order $\mathcal{O}(1/b_{k})$

[TABLE]

So for a higher intermediate value of $b_{k}$, which in the end also means more possibly pathways in the branching tree, it should be less favorable to change the current $\Delta b_{k}$ value. In order to create an excitation $\ket{m^{\prime}}$ with a probability proportional to the coupling coefficient $\lvert\bra{m^{\prime}}\hat{E}_{ij}\ket{m}\rvert$ this fact is included in the decision of the chosen branch and is achieved by the use of branch weights.

IV.5 Branch Weights

It is possible to take into account the “probabilistic weight” of each tree branches at a possible branching decision. As one can see in the left panel of Fig. 8, the starting $\Delta b_{i}=\pm 1$ branches each have one contribution of order $\mathcal{O}(1)$. For each branching possibility there is a resulting branch with opposite $\Delta b$ and weight of order $\mathcal{O}(b^{-1})$. However, it also depends on the end-segment determined by $d_{j}$, if a given branch can be chosen. The following branch weights

[TABLE]

with

[TABLE]

where $s_{k}(\pm 1)$ is the number of remaining switches (48), can be used to determine the probability of each $\Delta b$ branch to be chosen. The right panel of Fig. 8 shows the influence of the matrix element on the branching probabilities in the excitation range. By choosing the $\Delta b=-1$ with a probability $p_{-}=\frac{\zeta_{-}}{\zeta_{-}+\zeta_{+}}$ at the start of an excitation and in the excitation region choose to stay on the current $\Delta b$ branch according to

[TABLE]

the overall probability to choose the specific excitation $\ket{m^{\prime}}$ is given by

[TABLE]

With this choice of branching probabilities it is possible to retain an almost linear ratio between $p(m^{\prime}|m)$ and coupling coefficient amplitudes $\lvert\langle{m^{\prime}}|{\hat{E}_{ij}}|{m}\rangle\rvert$. Additionally, because of the $f(d_{j})$ and $g(d_{j})$ functions and inclusion of the remaining switches (48) in Eq. (50), this approach avoids dead-ends and thus choosing invalid excitations.

An important note on the matrix element calculation of single excitations: there are of course contractions of the two-body operator in Eq. (22), which contribute to the matrix element of a single excitation $\langle{m^{\prime}}|{\hat{H}}|{m}\rangle$. These contractions have to be taken into account in the “on-the-fly matrix element computation” and are explained in more detail in Appendix B or can be found in Ref. [(74)].

V Excitation Generation: Doubles

The generation of double excitation in the GUGA formalism is much more involved than single excitations and the detailed background on matrix element computation and weighted orbital choice can be found in Appendix C or Ref. [(74)] for conciseness of this manuscript. Here we will only present the general ideas involved in doubly excitation generation.

Depending on the ordering of the involved spatial orbitals of the one- and two-body generators, $\hat{E}_{ij}$ and $\hat{e}_{ij,kl}$, thirty different excitation types, involving different combinations of lowering (L) and raising (R) generators, can be identified and are listed in Table 6.

Some of them are equivalent, in the sense that they lead to the same excitations, such as the 7 single excitation (0a-0g) in Table 6, which reduce to the two distinct raising $\underline{R}\rightarrow\overline{R}$ and lowering $\underline{L}\rightarrow\overline{L}$ generators. The pictorial representation of these generators are shown in Fig. 9, where the ordering of orbitals is from bottom to top and arrows indicate the replacement of electrons.

The two-body operators $\hat{e}_{ij,kl}$, which contribute to single excitations, (0c-0g) in Table 6, are already accounted for in the single excitation matrix element calculation, see Sec. IV.1 and Appendix B.2. These also include the single overlap excitation (0b) and (0e) with two alike generator types.

Double excitation with a single overlapping index $j$ but two different generators (1a) and (1b) can be treated in a similar way to single excitations, with the same weighting functions (50) and classification of remaining switches (48), but with a change of generator type at the overlap site, $L\leftrightarrow R$. Double excitations with an empty overlap range $S_{1}$ (31) (3c0, 3d0, 3e0 and 3f0) can be calculated as the product of two single excitations (33). However, e.g. for excitation (3c0), the two-body operators $\hat{e}_{ij,kl}$ and $\hat{e}_{kj,il}$ contribute to the same Hamiltonian matrix element. We made the decision to treat these non-overlap excitations by using the corresponding two-body generators with a nonzero overlap range $S_{1}$, see App. C for more details.

For “proper” double excitations, we separate the excitation range $\min(i,j,k,l)\rightarrow\max(i,j,k,l)$ into the lower non-overlap range $S_{2}$ below the overlap range $S_{1}$ and the upper non-overlap range $S^{\prime}_{2}$ above $S_{1}$, as depicted in Fig. 5. We introduce the terminology of a full-start of mixed generators $\underline{R}\underline{L}$ and alike generators $\underline{R}\underline{R}/\underline{L}\underline{L}$, a semi-start corresponds to the segment types like $R\underline{R}$ or $R\underline{L}$, a semi-stop indicates generator combination like $L\overline{L}$ or $L\overline{R}$ and a full-stop is where both generators end on the same orbital, e.g. $\overline{L}\overline{L}$ or $\overline{R}\overline{L}$.

The excitation generation for doubles is again performed by choosing a valid path in a branching tree with modified rules in the overlap range $S_{1}$ of the double excitation. As can be seen by the restrictions for nonzero two-body matrix elements, Eq. (35-43), the allowed $\Delta b$ values in $S_{1}$ are now $\pm 2$ and [math]. This leads to new elements of the branching tree in $S_{1}$, which are shown by the example of alike raising and mixed generators in Fig. 10, where vertical lines indicate the new $\Delta b=0$ branch, and left (right) going lines in $S_{1}$ correspond to $\Delta b=\pm 2$. The rules for the intermediate elements $RR,LL$ and $RL$ are the same for all combinations of generators.

The calculation of the remaining switch possibilities (48) essentially is the same as for single excitations, except they are calculated for each segment, $S_{2},S_{1}$ and $S^{\prime}_{2}$ of the excitation separately. In $S_{1}$ a $\Delta b_{k}=-2$ branch can switch at $d_{k}=1$, a $\Delta b_{k}=+2$ at $d_{k}=2$ and the $\Delta b_{k}=0$ branch at both open-shell step-values

[TABLE]

The remaining switches in $S_{2}$ are calculated up until the index of the start of $S_{1}$, as, similar to single end segments, e.g. $\overline{R}$, there are the restrictions for nonzero matrix elements for semi-start segments, e.g. $R\underline{L}$, to guarantee the total spin is conserved. Similarly, for the end of the overlap range, depending on the step-value at e.g. $L\overline{L}$, the mentioned restrictions apply so the remaining switches (54) are calculated until the start of $S^{\prime}_{2}$.

To relate $p(m^{\prime}|m)$ to the generator matrix element $\langle{m^{\prime}}|{\hat{e}_{ij,kl}}|{m}\rangle$ we again use branching weights to determine which paths of the tree are chosen. For a full-start $\underline{RL}(i)$ into full-stop $\overline{RL}(j)$ excitation the weights of the different $\Delta b$ branches in terms of the intermediate $b$-values and remaining switch possibilities are

[TABLE]

with $f(d_{j})$ and $g(d_{j})$ given by Eq. (51). We bias towards the $\Delta b_{k}=0$ branch at the start of the excitation range with

[TABLE]

depending if $d_{i}=\{1,2\}$ and weight to stay on the current $\Delta b_{k}$ excitation branch in $S_{1}$ with

[TABLE]

For a full-start into semi-stop excitation, e.g. $\underline{RL}\rightarrow\overline{L}R\rightarrow\overline{R}$, the weights of the branches in the overlap region $S_{1}$ are given by

[TABLE]

where $\zeta_{\pm 1}(j)$ are the single weights (50) for the non-overlap region $S^{\prime}_{2}$ at the end of the excitation, evaluated with the $b_{j}$ and $s_{j}(\pm 2)$ values at the semi-stop. The biasing function towards a certain branch at the beginning of an excitation and to stay at a chosen $\Delta b$ branch are the same as (58) and (59) and in the non-overlap region $S^{\prime}_{2}$ the single excitation weights and biasing factors (50, 52) apply.

The weighting functions in the non-overlap region $S_{2}$ for a semi-start into full-stop excitation, e.g. $\underline{R}(i)\rightarrow\underline{L}R(j)\rightarrow\overline{R}\overline{L}(k)$, are given by

[TABLE]

with $\Sigma_{x}$ being the weights of the full-stop excitation (55) evaluated with the $b_{j}$ and $s_{j}(\pm 2)$ values at the start of the overlap region $j$. The biasing function for the start and staying probabilities are the same as in the single excitation case (52) evaluated with $\sigma_{\pm 1}$ instead of $\zeta_{\pm 1}$.

For a “full” double excitation, e.g. $\underline{R}(i)\rightarrow\underline{R}R(j)\rightarrow\overline{R}R(k)\rightarrow\overline{R}(l)$, the weights and biasing functions for the first non-overlap $S_{2}$ region $i\rightarrow j-1$ are the same as for the semi-start into full-stop excitation $(\ref{eq:semistart-fullstop-weights})$, but evaluated with the full-start into semi-stop weights $\Sigma_{x}$(60). In the overlap region $S_{1}$, $j\rightarrow k-1$, the weights and biasing functions are the same as for full-start into semi-stop excitations (60), where the $f(d_{k}),g(d_{k})$ and $\zeta_{\pm 1}$ functions are evaluated at the semi-stop index $k$ now. And finally for the final non-overlap region $S^{\prime}_{2}$, $k\rightarrow l-1$, the weights and biasing functions for single excitations (50, 52) apply.

By using this biasing we ensure to create a valid spin conserving excitation, avoid ending up in a dead-end of the branching tree and create excitations with a probability $p(m^{\prime}|m)$ proportional to the coupling coefficient magnitude $\lvert\langle{m^{\prime}}|{\hat{e}_{ij,kl}}|{m}\rangle\rvert$. The used weight functions are set up before an excitation in terms of the $b_{k}$ and remaining switch possibilities with, if necessary, the precomputed switch possibilities for the remaining overlap and non-overlap contributions and $\Delta b$ conditions. It is not necessary to recompute the whole setup at each step of the excitation. The computational effort to set up this weight objects, as it needs the information of the remaining switches, is $\mathcal{O}(n)$, in the worst case of an excitation spanning the whole orbital range. An analysis of the increase in computational effort of the GUGA-FCIQMC method compared to the SD based implementation can be found below.

VI Histogram based Timestep Optimization

Due to the increased connectivity of CSFs compared to SDs, the generation probability, $p(m^{\prime}|m)$, to spawn a new walker on state $\ket{m^{\prime}}$ from an occupied CSF $\ket{m}$, is in general much lower than between SDs. An efficient sampling of the off-diagonal Hamiltonian matrix elements and stable dynamics of a simulation, demand the quantity $\Delta\tau\lvert H_{m^{\prime}m}\rvert/p(m^{\prime}|m)$ to be close to unity. In the original determinant-based FCIQMC algorithm this is ensured by a dynamically adapted timestep $\Delta\tau(t)$, taking on the value of the “worst-case” $p(m^{\prime}|m)/\lvert H_{m^{\prime}m}\rvert$ ratio encountered during a simulation.

However, due to the large number of possible connections between CSFs, this causes the timestep to drop dramatically. At the same time a tiny spawning probability $p(m^{\prime}|m)$ means that these problematic excitation only happen a minuscule fraction of times compared to more “well-behaved” excitations. Through the timestep, the global dynamics of all the walkers are affected by possibly only one ill-sampled excitation with a large $\lvert H_{m^{\prime}m}\rvert/p(m^{\prime}|m)$ ratio. The optimized excitation generation mentioned in the sections above, ameliorates this issue, but still cannot avoid the inherent “connectivity problem” of a CSF based implementation. If we store all $\lvert H_{m^{\prime}m}\rvert/p(m^{\prime}|m)$ of all successful excitation attempts in a histogram of certain bin width, we can see that the majority of excitations are well represented by the optimized generation probability, see Fig. 12.

The SD based method has a fast exponential decaying tail. This is the reason the “worst-case” timestep adaptation does not cause any problems for the original FCIQMC implementation. The GUGA implementation on the other hand, especially in the unoptimized version (uniform choice of branching possibilities and now weighting according to the molecular integrals), has a very slow decay and much larger maximum $\lvert H_{m^{\prime}m}\rvert/p(m^{\prime}|m)$ ratios, over 10000 in the N2 example shown in Fig. 12 (not displayed for clarity). The optimized CSF excitation scheme, explained above, greatly improves the $p(m^{\prime}|m)$ to $\lvert H_{m^{\prime}m}\rvert$ relation, but expectedly behaves worse than the SD based method. The timestep obtained with the “worst-case” optimization are given in Table 7.

To avoid this hampering of the global dynamics by a few ill-behaved excitations, we implemented a new automated timestep adaptation by storing the $\lvert H_{m^{\prime}m}\rvert/p(m^{\prime}|m)$ ratios off all successful excitation attempts in a histogram, and setting the timestep $\Delta\tau$ to ensure $\Delta\tau\lvert H_{m^{\prime}m}\rvert/p(m^{\prime}|m)\leq 1$ for a certain percentage of all excitations. The results of this “histogram-tau-search” are listed in Table 7 for a SD based and unoptimized (vanilla) and optimized GUGA-FCIQMC implementation for simulations of the nitrogen dimer at equilibrium geometry in a cc-pVDZ basis. For an SD-based implementation there is not much difference between the two approaches. Similar, for the vanilla GUGA implementation, due to the slow decaying tail in the histograms, see Fig. 12. There is a two order of magnitude difference between the SD-based and the unoptimized GUGA-based timestep, which in practice would make the GUGA-FCIQMC implementation useless. However, with the optimized CSF excitation generation, the histogram-based $\Delta\tau_{h}$-adaptation yields a timestep two orders of magnitude larger than the “worst-case” $\Delta\tau_{w}$-optimization. The obtained $\Delta\tau_{h}$ is still half that of the SD-based FCIQMC, but due to a smaller Hilbert space size, and possibly faster convergence for spin-degenerate systems, this makes the GUGA-FCIQMC applicable for real systems.

VII Results and Discussion

VII.1 Nitrogen Atom

To benchmark the GUGA-FCIQMC implementation we first investigated the nitrogen atom. The ground state configuration of N is 1s22s22p3 with the 3 electrons in the p-shell forming a $S=3/2$ quartet ${}^{4}S^{o}$ state. The first excited state is the $S=1/2$ ${}^{2}D^{o}$ doublet, $2.384$ eV above the ground state Kramida et al. (2018); Gallagher and Moore (1993), with spin-orbit effects neglected. This setup of a half-integer high-spin ground state with low-spin excited state is the prime playground of the GUGA-FCIQMC method. Previous spin-adapted implementations in FCIQMC, using half-projected and projected Hartree-Fock (HPHF) states Smeyers and Doreste-Suarez (1973), are only applicable to an even number of electrons. At the same time, restricting the total $m_{s}$ quantum number to target an excited state only works if the low-spin state is the ground state with excited states being high-spin, since the high-spin ground state also contains contributions of energetically lower $m_{s}$ states, causing the projective FCIQMC to converge to the latter one.

We prepared all-electron ab-initio Hamiltonians with MOLPRO Werner et al. (2012, 2015) for N in a cc-pV$n$Z basis set, with $n$ = D, T, Q, 5 and 6. The maximal symmetry point group in MOLPRO is D2h and thus the much larger SO(3) symmetry of N gets reduced to the one-dimensional irreps of D2h. The $S=3/2$ quartet with singly occupied 2p orbitals belongs to the irrep Au. While the $S=1/2$ doublet splits into one Au state with three open-shell 2p orbitals and three states belonging to B1u, B2u and B3u with one doubly occupied and one open-shell 2p orbital, see Fig. 13.

We calculated the quartet-doublet, ${}^{4}S^{o}-{{}^{2}D^{o}}$, spin-gap with the spin-adapted i-FCIQMC method (GUGA-FCIQMC) for basis sets up to cc-pV6Z and compared our results to unrestricted coupled cluster singles and doubles with perturbative triples ((U)CCSD(T)) and FCI calculations up to cc-pVTZ obtained with MOLPRO Knowles and Handy (1984, 1989); Knowles, Hampel, and Werner (1993, 2000) and experimental results Kramida et al. (2018); Gallagher and Moore (1993). The CCSD(T) calculations are based on restricted open-shell Hartree-Fock Roothaan (1960) (ROHF) orbitals, which for the $S=1/2$ state are only possible to be done for the Biu, $i=1,2,3$, states. Although GUGA-FCIQMC calculations for the Biu irreps yield the same results as for the $S=1/2$ Au state, the CCSD(T) results are far off the FCI results and the experimental gap, due to the multi-reference character of these states.

The results are given in Table 8 with a complete basis set (CBS) extrapolation given by a two-parameter inverse cube fit Helgaker et al. (1997)

[TABLE]

using $n$ = T, Q and 5. The GUGA-FCIQMC CBS results shows excellent agreement with the experimental value within chemical accuracy, while the CCSD(T) calculations are not able to obtain the correct result, due to the multiconfigurational character of the ${}^{2}D^{o}$ excited state.

We also calculated the ionization potential (IP) of the nitrogen atom in the CBS limit with GUGA-FCIQMC and compared our results to CCSD(T) calculations and experimental data. The ground state of the N+ cation is the $S=1$ triplet ${}^{3}P_{0}$ state. The results from GUGA-FCIQMC and CCSD(T) calculations up to cc-pV6Z basis set are shown in Table 8. Since CCSD(T) Knowles and Handy (1984, 1989); Knowles, Hampel, and Werner (1993, 2000) can treat both the ${}^{4}S^{o}$ and ${}^{3}P_{0}$ well, coupled cluster results and GUGA-FCIQMC CBS limit values, using $n$ = Q, 5 and 6, agree within chemical accuracy with experimental values.

VII.2 Nitrogen Dimer

The breaking of the strong triple bond of N2 is accompanied by a change of single-reference to multiconfigurational character of the electronic structure and the concomitant strong electron correlation effects pose a difficult problem for quantum chemical methods. The ground state of the nitrogen molecule at equilibrium bond distance, $r_{0}\approx 2.1a_{0}$, is a singlet ${}^{1}\Sigma_{g}^{+}$, where all bonding molecular orbitals (MOs) formed from the 2p atomic orbitals (AOs) of the constituent N atoms are doubly occupied

At large bond distances the ground states of the $S=0,1,2$ and $3$ states are degenerate, since the coupling of the independent nitrogen atoms $A$ and $B$, ${}^{4}S^{o}_{A}\otimes^{4}S^{o}_{B}$, are all degenerate. We calculated the dissociation energy of N2 as the difference of the ${}^{1}\Sigma_{g}^{+}$ N2 ground state at equilibrium geometry $r_{0}=2.074a_{0}$ and the ${}^{7}\Sigma_{u}^{+}$ state at $r=30a_{0}$ in the cc-pV$n$Z basis set, up to $n=5$, with four core electrons frozen and performed a CBS limit extrapolation using Eq. (65) with the $n$ = T, Q and 5 results. The results are shown in Table 9 with CCSD(T) results obtained with MOLPRO Werner et al. (2012, 2015); Knowles, Hampel, and Werner (1993, 2000) and compared with experimental results Huber and Herzberg (2013), which are corrected to remove scalar relativistic, spin-orbit and core correlation effects according to Refs. (87; 88). We also checked the convergence of the $r=30a_{0}$ results with the independent atom calculations with a frozen core and found excellent agreement. Both the GUGA-FCIQMC and CCSD(T) results agree with experimental values within chemical accuracy.

To investigate the correct accounting of core and core-valence correlation effects, we performed all-electron GUGA-FCIQMC and CCSD(T) calculations in the cc-pCV$n$Z basis set and calculated the dissociation energy of N2 as the difference of the independent nitrogen atom 4So ground state results in the same basis set and the ${}^{1}\Sigma_{g}^{+}$ N2 ground state at equilibrium geometry $r_{0}=2.074a_{0}$, as $E_{diss}=2E_{atom}-E_{dimer}$. The results are shown in Table 9 and the CBS limit extrapolations of the N2 dissociation energy agree within chemical accuracy with experimentHuber and Herzberg (2013) for both the GUGA-FCIQMC and CCSD(T) calculations. We also performed a counterpoise correction Boys and Bernardi (1970), but found the basis set superposition error to be negligibly small.

To show the improved convergence behavior of the spin-adapted FCIQMC method for systems with near-degenerate spin-eigenstates, we calculated the gap of the singlet ${}^{1}\Sigma_{g}^{+}$ ground-state to the triplet ${}^{3}\Sigma_{u}^{+}$, quintet ${}^{5}\Sigma_{g}^{+}$ and septet ${}^{7}\Sigma_{u}^{+}$ excited states of N2—which are all degenerate at dissociation—for the equilibrium bond distance $r=2.118\,a_{0}$111In the cc-pVDZ basis set. and two stretched geometries $r=4.2\,a_{0}$ and $r=6.0\,a_{0}$ in a cc-pVDZ basis set with the original SD-based and GUGA-FCIQMC method. Figure 14 shows the gaps between the ground state and three excited states as a function of the total walker number compared with DMRG reference results Chan and Head-Gordon (2002); Sharma and Chan (2012). Since the energy of the spin states are ordered according to their total spin quantum number, it is possible to obtain the spin-gaps in the original determinant based FCIQMC method by restricting the $m_{s}$ quantum number alone. At cc-pVDZ equilibrium bond distance $r=2.118\,a_{0}$ the determinant based and spin-adapted FCIQMC implementation are equivalent in their convergence behavior of the spin-gaps w.r.t. the walker number. However, as the bond distance increases, and thus the spin-gaps decrease, the spin-adapted FCIQMC implementation shows a dramatically improved convergence, especially for the singlet-quintet and singlet-septet gaps, where around an order of magnitude fewer walkers are necessary to obtain the same accuracy as the SD-based FCIQMC method.

The second common option to obtain spin-gaps with the FCIQMC method is based on HPHF functions, which allow targeting spin-states with an even (singlet, quintet, …) or odd (triplet, septet, …) total spin $S$, allowing to obtain the singlet-triplet gap in the case of N2. Figure 15 shows the relative error of the singlet triplet gap, obtained with the HPHF based and spin-adapted FCIQMC implementation with $N_{w}^{tot}=10^{7}$ as a function of the singlet-quintet gap, compared to DMRG reference results Chan and Head-Gordon (2002); Sharma and Chan (2012) on a double logarithmic scale. As both even-spin singlet and quintet states belong to the same spatial point group irrep Ag, the HPHF solution is spin-contaminated by an increasing amount for a decreasing singlet-quintet gap. This fact prohibits the HPHF-based FCIQMC implementation to obtain the correct singlet-triplet gaps for increasing bond distance for the nitrogen dimer.

VII.3 Computational Effort and Scaling of GUGA-FCIQMC

To analyze the additional computational cost associated with the GUGA-based CSF implementation in FCIQMC, we compare the time per iteration, $t_{iter}$, and timestep, $\Delta\tau$, with the original SD-based FCIQMC method for the nitrogen atom and dimer, mentioned above. Since FCIQMC is formally linear-scaling with the walker number $N_{w}$Booth, Smart, and Alavi (2014) we removed the bias of walker number differences by comparing the time per iteration and per walker.

The left panel of Fig. 16 shows the timestep $\Delta\tau$ obtained with the histogram based optimization, see Sec. VI, for N2 at $r=4.2\,a_{0}$ vs. the cardinal number $n$ of the cc-pV$n$Z basis set. As expected the usable timestep in the SD-based simulation is higher compared to the CSF-based calculation, with roughly twice the possible $\Delta\tau$. However, rather surprisingly the difference between the two decreases with increasing basis set size. The right panel of Fig. 16 shows the time per iteration and walker for the same simulations. The additional computational cost of the GUGA implementation roughly doubles the time per iteration compared to the original FCIQMC method. While there seems to be a steeper increase with increasing basis set size for the CSF-based implementation, it is nowhere near the formally $\mathcal{O}(n)$ cost, with $n$ being the number of orbitals, mentioned above. In total, with twice the timestep and twice the time per iteration, the spin-pure GUGA implementation amounts to a fourfold increase in computational cost compared to the original SD-based FCIQMC method. 222All simulations for this comparison were performed on identical 20 core Intel Xeon E5-2680 nodes with 2.8GHz clock rate, 20MB cache and 128GB memory.

To examine the scaling in more detail, a least-squares fit to the polynomial $f(n)=a+b\cdot n^{c}$, with 3 parameters $a$, $b$ and $c$, was performed on the available data points with $n$ being the cardinal number of the basis set. The lines in Fig. 16 represent this fit for the timestep $\Delta\tau(n)$ and time per iteration $t_{iter}(n)$, as a function of the cardinal number $n$ of the basis set. The results of the least-squares fit for the determinant- and CSF-based calculations are shown in Table 10.

The scaling of the decrease in the possible timestep $\Delta\tau$ is almost less than a factor of $n$ smaller in the CSF based implementation and the increase of the time per iteration $t_{iter}$ less than $n$ larger compared to the determinant based implementation. However, the combination of these two effects causes the spin-adapted FCIQMC implementation to scale by an additional factor of $\approx\mathcal{O}(n^{1.3})$ for this specific system, compared to the original SD-based FCIQMC method.

Table 11 shows the averaged timestep and time per iteration ratios between GUGA- and SD-based simulations for the nitrogen atom and dimer. Compared to the CSF-based FCIQMC the maximum possible timestep in the original determinant based implementation is larger by a factor of $2.55$ to $3.87$ and the time per iteration is smaller by a factor of $0.63$ to $0.90$. The combination of these effects result in a slow down by a factor of $2.8$ to $5.0$ of the spin-adapted FCIQMC method.

VII.4 The Cobalt Atom

As with most open-shell transition metals, the cobalt atom has a high-spin ground state, due to Hund’s first rule. This prohibits the calculation of the spin-gap to low-spin excited states by restriction of the $m_{s}$ quantum number, as inevitably these excited state calculations will converge to the high-spin ground state in the projective procedure of FCIQMC.

We compare our results to coupled cluster calculations, which are not so easily applicable, due to the multireference character of the excited states of these systems.

The ground state electronic configuration of the neutral cobalt atom is [Ar]$3s^{2}3p^{6}3d^{7}4s^{2}$ and is a quartet ${}^{4}F$ state. We calculated the spin gap to the first doublet excited state ${}^{2}F$ with the [Ar]$3s^{2}3p^{6}3d^{8}4s$ configuration with the GUGA-FCIQMC method, correlating 17 electrons in all available orbitals. We employed an ANO basis set Almlöf and Taylor (1991) with primitive contractions corresponding to a comparable V$n$ZP basis with $n$ = D, T and Q and the full and completely uncontracted primitive ANO basis set with 2nd order Douglas-Kroll scalar relativistic corrections Reiher (2006). The ANO molecular integral files were computed with MOLCAS Aquilante et al. (2016). We also prepared ab-initio integrals with an augmented correlation consistent core-valence basis set with 2nd order Douglas-Kroll scalar relativistic corrections Reiher (2006), aug-cc-pwCV$n$Z-DK (denoted as cc-basis in Table 12), up to $n$ = Q. The cc-basis molecular integrals were computed with MOLPRO Werner et al. (2012, 2015). We also performed 2nd order complete active space perturbation theory Malmqvist, Rendell, and Roos (1990); Andersson, Malmqvist, and Roos (1992) (CASPT2) calculation on the ANO basis set with MOLCAS and CCSD(T) calculation in the cc-basis with MOLPRO.

The starting orbitals for the Co ${}^{2}F$ and ${}^{4}F$ calculation with FCIQMC were CASSCF Olsen (2011); Hegarty and Robb (1979) orbitals with the $1s^{2}2s^{2}2p^{6}3s^{2}3p^{6}$ orbitals frozen, 9 active electrons in the active space of $4s,3d,4p,5s$ and $4d$, CAS(9,15) and further orbitals being virtuals, see Fig. 17. The CASSCF calculations were performed with MOLCAS Aquilante et al. (2016) and MOLPRO Werner and Knowles (1985); Knowles and Werner (1985). Similar to the nitrogen atom the SO(3) symmetry of Co is reduced to the D2h symmetry implemented in MOLPRO and MOLCAS. We chose the B1g irrep for the ${}^{2}F$- and the Ag irrep for the ${}^{4}F$-state.

Similar to the nitrogen atom the odd number of electrons and high-spin ground state to low-spin excited state setup makes previous spin-pure methods implemented in FCIQMC not applicable. However, as the results in Table 12 show, the GUGA-FCIQMC implementation is able to provide energies within chemical accuracy close to the experimental result Kramida et al. (2018); Sugar and Corliss (1985); Pickering and Thorne (1996). For both GUGA-FCIQMC and CASPT2, the CBS limit extrapolation of the spin-gap, using the VTZP and VQZP results for Eq. (65), in the ANO basis set agree within 1 kcal/mol (chemical accuracy) with the experimental result.

For the aug-cc-pwCV$n$Z-DK we performed separate CBS limit extrapolations of the 2F and 4F ground state energy, using the Hartree-Fock energy of a $n=5$ calculation and a two-point extrapolation of the correlation energy, according to Eq. (65), using the $n=$T and Q data points. The resulting estimated spin-gap lies approximately $0.0024\,E_{h}\approx 1.519$ kcal/mol above the experimental result, see Table 12.

Similar to the spin gap of nitrogen, see Sec. VII.1, coupled cluster is not able to provide correct results of the doublet ${}^{2}F$ state of cobalt. The CCSD(T) calculations are based on ROHF orbitals and the valence electronic configuration of the ${}^{2}F$ state, $3d^{8}4s$, enforces the $4s$ orbital to be singly occupied with all the $d$-electrons being in a closed shell conformation. This obviously violates Hund’s rule and thus the CCSD(T) results give a too high energy for the ${}^{2}F$ state and thus the spin-gap is immensely overestimated, see Table 12. Further investigations are being conducted on the performance of CCSD(T) by varying the reference determinant, and will be reported elsewhere.

VIII Discussion

The spin-adapted FCIQMC implementation was tested and benchmarked for the nitrogen atom and dimer, where excellent agreement with exact results—where available—and other quantum chemical methods was observed. We found that the additional computational cost associated with the more complicated and highly connected Hilbert space of CSFs is manageable and applications of this approach for large basis sets was demonstrated; eradicating the severe limitations of previous spin-adapted approaches in general and in FCIQMC in particular.

The validity of the approach was proven and the direct targeting of specific spin states is possible; enabling us to obtain results previously not accessible to the FCIQMC method. These are gaps of high-spin ground and low-spin excited state systems with an odd number of electrons and the excitation energies within an explicit spin symmetry sector. For system with near-degenerate spin-eigenstates we observe an accelerated convergence of spin-gap results with respect to the total walker number in the spin-adapted FCIQMC implementation.

However, the additional scaling with the number of spatial orbitals in the GUGA-FCIQMC method, starts to become relevant for large basis set expansions, limiting the applicability, where the SD based implementation remains preferable. The increased connectivity of a spin-pure basis reduces the generation probabilities in the spawning step of the FCIQMC method and thus limits the possible timestep of a simulation and causing stability issues in the sampling process. In this regard, the scope of application of this method is to target specific, interesting spin states, which allows a clearer chemical and physical interpretation of results. As a consequence, more insight in chemical processes governed by the interplay of different spin states is possible.

IX Conclusion and Outlook

The efficient usage of a spin-adapted basis in FCIQMC has been made possible within the (graphical) unitary group approach (GUGA) and the severe limitations of previous implementations have been overcome. When formulated in such a basis, simulations conserve the total spin quantum number and the Hilbert space size of the problem is reduced. As another positive consequence, targeting specific many-body subspaces of the Hamiltonian and getting access to their excitation energies is possible, and thus one is able to study phenomena governed by the interplay of different—even degenerate—spin sectors. Additionally, the use of a spin-adapted basis improves the convergence of the projective FCIQMC method, for systems with near-degenerate spin states.

We benchmarked the spin-adapted FCIQMC method and compare results with other computational approaches, for the nitrogen atom and dimer, where we find excellent agreement with reference results, when available. For the nitrogen atom we obtain the spin gap of the 4So ground- and 2Do excited state and the ionization potential, and the dissociation energy of the nitrogen dimer within chemical accuracy to experiment. We apply the method to study the 3d-transition metal cobalt, targeting properties, which defy a simple single-reference description. For cobalt, the spin-gap of the high-spin ground state (single-reference wavefunction) and low-spin excited state (multi-reference wavefunction) was determined within chemical accuracy to experiment.

This spin-adapted implementation brings FCIQMC en par with many other quantum chemical methods, which already utilize the inherent total spin conservation of nonrelativistic, spin-independent molecular Hamiltonians.

To combine the spin-adapted FCIQMC with the stochastic CASSCF method, the final missing piece is the spin-pure implementation of an efficient sampling of reduced density matrices Overy et al. (2014), which would enable us to solve active spaces of unprecedented size in a spin-pure fashion, extending even further the applicability of the method. Unfortunately, the sampling of RDMs in the spin-adapted formulation based on the GUGA is unfortunately a highly non-trivial task. Although there is no theoretical problem of density matrices in the unitary group formalism Gould, Paldus, and Chandler (1990); Paldus and Gould (1993); Shepard (2006); Polyak, Bearpark, and Robb (2018), from a practical standpoint there is. Due to the increased connectivity within a CSF basis and the possibility of generators with different spatial indices contributing to the same density matrix element, there is a large overhead involved in sampling RDMs in the spin-adapted FCIQMC method. However, we are optimistic to solve these problems in due time, which would allow us to use GUGA-FCIQMC as a spin-pure FCI solver in the stochastic CASSCF method Li Manni, Smart, and Alavi (2016). This would enable us to solve active spaces of unprecedented size in a spin-pure fashion, extending even further the applicability of the method. Furthermore, the unitary group formalism is extendable to spin-dependent operators Drake and Schlesinger (1977); Kent, Schlesinger, and Drake (1981); Gould and Chandler (1984); Kent and Schlesinger (1990); Gould and Paldus (1990); Gould and Battle (1993); Kent and Schlesinger (1994); Li and Paldus (2014); Yabushita (2014), and an extension of the spin-adapted FCIQMC method to this approach is currently investigated to enable us to study systems with spin-orbit coupling and explicit spin dependence. Along this line, another interesting problem to be investigated is the application of GUGA-FCIQMC to the two-dimensional $t$-$J$ and Heisenberg models.

Appendix A CSF Excitation Identification

Efficiently identifying the difference between two given CSFs and the type of excitation (generator types), listed in Table 6 of the main text, is crucial for an optimized matrix element calculation. For CSFs this operation is more involved compared to Slater determinants. This is because not only occupancy differences but also changes in the singly occupied orbitals (different spin-couplings) must be taken into account, as they can also lead to non-zero coupling coefficients. The defining difference for the excitation is the difference in spatial occupation numbers. The step-values, $d_{i}=\{0,1,2,3\}$, of the spatial orbitals of a CSF are efficiently encoded by two bits per spatial orbital

[TABLE]

in an integer of length $2n$. This is equivalent to the memory requirement of storing the occupied spin-orbitals of a Slater determinant. The spatial occupation difference, $\lvert\Delta n\rvert$, can be efficiently obtained by shifting all negatively spin-coupled, $d_{i}=2:10$, to the right, and computing the bit-wise xor-operation on two given CSFs and counting the number of set bits in $\lvert\Delta n\rvert$, e.g. by the Fortran 2008 intrinsic popcnt:

[TABLE]

With $\Sigma\lvert\Delta n\rvert$ we can identify the excitation level, which would be a single excitation from orbital $4$ to $1$ in the example above, and $\Delta n$ gives us information, in which spatial orbitals electrons got removed or added. For CSFs the orbital occupation difference alone is not enough to completely identify an excitation between two CSFs, since for excitations of exchange type, involving $\underline{R}\underline{L}$ and $\overline{R}\overline{L}$ generators, there can be a change in the spin-coupling, without an actual change in orbital occupation. So additionally we also need information of the step-vector difference, $\Delta d$, which is just obtained by the xor-operation on the bit-representation of two given CSFs:

[TABLE]

In this example it can be seen, that the $\lvert\Delta n\rvert$ information alone would lead us to believe a single excitation connects $m$ and $m^{\prime}$, but this is not compatible with the change in step-vector at orbital $2$. So in addition, we need to determine if there are step-vector changes below the first, $\Delta d_{b}$, or above the last, $\Delta d_{a}$, occupation change in $\Delta n$. This can be done efficiently with the Fortran 2008 intrinsic bit-operations, leadz(I) (trailz(I)), which give the number of leading(trailing) zeros in integer I. The case that there are only step-vector changes, $\Delta d$, within, the first and last $\Delta\neq 0$ cases, is encoded by $\Delta d_{b}=\Delta d_{a}=0$.

$\Sigma\lvert\Delta\rvert n>4$ indicates a higher excitation than double, so the two CSFs are not possibly connected by a single Hamiltonian application and can be disregarded. The non-zero Hamiltonian matrix elements can be identified by following combinations of $\Delta n$ and $\Delta d$:

$\Sigma\lvert\Delta n\rvert=0\;\&\;\Delta d\neq 0$:

This combination indicates, that there is no difference in the occupation number between two CSFs $m$ and $m^{\prime}$, but a change in the spin-coupling of the singly occupied orbitals. Only a mixed generator $\underline{R}\underline{L}\rightarrow\overline{R}\overline{L}$ generator combination, corresponding to the type (2c) in Table 6, can lead to those types of excitations. Details on the matrix element calculation in general are discussed below.

$\Sigma\lvert\Delta n\rvert=2\;\&\;\Delta d_{b}=\Delta d_{a}=0$:

This combination indicates a regular single excitation and the order of the removed and added electron determines the type of generator $\hat{E}_{ij}$, corresponding to type (0a) in Table 6,

[TABLE]

$\Sigma\lvert\Delta n\rvert=2\;\&\;\Delta d_{b}\neq 0\;\text{or}\;d_{a}\neq 0$:

This indicates step-vector changes above or below the occupation differences, which identifies a mixed start $\underline{R}\underline{L}$ or end $\overline{R}\overline{L}$ segment. Again the order of the orbital occupation and step-vector changes below or above $\Delta n$ identifies the type of excitation

[TABLE]

with $\Delta n_{\pm 1}=\Delta n=\pm 1$ and the reference to the entries of Table 6.

$\Sigma\lvert\Delta n\rvert=4\;\&\Delta d_{b}=\Delta d_{a}=0$:

In this case it is necessary to have $\Delta d_{b}=\Delta d_{a}=0$, otherwise this would indicate more than a double excitation, which would lead to a vanishing Hamiltonian matrix element. Again the order of the occupation differences gives information on the type of generators involved. The following combinations are identifiable only with $\Delta n$ (with reference to the entries of Table 6)

[TABLE]

where, e.g. $\underline{L}L/L\underline{L}$, indicates that the order of indices of equivalent generators, $\hat{e}_{il,jk}/\hat{e}_{ik,jl}$ see Fig. 18, can not be determined by $\Delta n$ alone. In the case of alike generators $RR(LL)$ this order does have influence on the sign of the matrix element Shavitt (1981), see below.

There are combinations of occupation differences where additionally the step vector differences have to be checked, since there are multiple two-body operators $\hat{e}_{ij,kl}$ possible, which can lead to the same excitation:

[TABLE]

where every second case is only possible if there are no step-vector differences, $\Delta d=0$, between the second and third occupation difference and reference to Table 6. The equivalence of these generator combinations can be seen in Fig. 19. However, even with no $\Delta d$ difference between the second and third occupation difference both generator combinations still contribute to the Hamiltonian matrix element.

With this method the distinct excitation types, listed in Table 6 of the main text, can be efficiently identified. We would like to mention that this type of excitation identification does not take into account the cumulative $b_{k}$ restrictions, $\lvert b_{k}-b^{\prime}_{k}\rvert\leq 1$ for single- and $\lvert b_{k}-b^{\prime}_{k}\rvert\leq 2$ for double-excitations, for non-zero matrix elements, which currently has to be directly accounted for in the subsequent matrix element calculation. However, we plan to make use of this property in future optimizations of our implementation, as it severely limits the number of nonzero coupling coefficients.

Appendix B Detailed Matrix Element Evaluation in the GUGA

In this section we explain the efficient matrix element calculation in a spin-pure CSF basis based on the GUGA approach in more detail.

B.1 Diagonal matrix elements

The diagonal matrix element for a GT state, $\ket{m}$ of the spin-free Hamiltonian, given in the main text, is given by the sum of the one-body matrix elements, $\langle{m}|{\hat{H}_{0}}|{m}\rangle$, with $\hat{H}_{0}=\sum_{ij}t_{ij}\,\hat{E}_{ij}$, and the two-body contribution, $\langle{m}|{\hat{H}_{1}}|{m}\rangle$, with $\hat{H}_{1}=\sum_{ijkl}V_{ijkl}\,\hat{e}_{ij,kl}$. The matrix elements of the weight generators, $\hat{E}_{ii}$, are simply just the occupancy of orbital $i$ in state $\ket{m}$

[TABLE]

So the one-body contribution is given by

[TABLE]

where $d_{i}$ is step-value of spatial orbital $i$ in $\ket{m}$. The two-body contributions are a bit more involved. Let’s consider the different cases:

$i=j=k=l$: reduces to the sum of doubly occupied orbitals

[TABLE]

since $\hat{E}_{ii}^{2}=\hat{E}_{ii}$ for $d_{i}={1,2}$.

$i=j\neq k=l$: only the $\hat{E}_{ii}\hat{E}_{jj}$ part of $\hat{e}_{ij,kl}=\hat{E}_{ij}\hat{E}_{kl}-\delta_{jk}\hat{E}_{il}$ remains, which reduces to a product of occupation numbers

[TABLE]

The last relation comes from the fact that the sums are invariant under $i$, $j$ exchange.

$i=l\neq j=k$: Also the exchange integral terms $V_{ijji}$ contribute to the diagonal matrix elements, if the excitation $\bra{m}\hat{E}_{ij}\hat{E}_{ji}\ket{m}$ leads to the same CSF. The calculation of these matrix elements depends on the step-values between $i$ and $j$ and are obtained by Shavitt’s graph rules Shavitt (1981). The matrix elements between two CSFs for a double excitation are given by the product

[TABLE]

with $S_{2}=(i,j)\cup(k,l)-S_{2}$ being the non-overlap range and $S_{1}=(i,j)\cap(k,l)$ being the overlap region of the indices of the involved generator $\hat{e}_{ij,kl}$. The one-body segment values $W(Q_{p};d_{p},d^{\prime}_{p},\Delta b_{p},b_{p})$ can be found in Table 3 in the main text and the two-body segment values $W_{x}(Q_{p};d_{p},d^{\prime}_{p},\Delta b_{p},b_{p})$ in Ref. [(69)]. In the case relevant for diagonal terms the matrix elements depending on the beginning $d_{i}$ and end $d_{j}$ step-values are given in Table 13 with

[TABLE]

Unfortunately this requires the consideration of all step-vector and b-values between $i$ and $j$ to calculate the diagonal matrix element.

[TABLE]

The first term $\frac{1}{2}n(d_{i})n(d_{j})$ accounts for the first singlet coupled $x=0$ matrix elements in the Table 13 and $X(i,j)$ accounts for the triplet coupled $x=1$ matrix elements. And only yields a contribution if both $d_{i}$ and $d_{j}$ are either 1 or 2

[TABLE]

and the rest gets accounted by the product of occupation numbers $n(d_{i})n(d_{j})$. In total the diagonal Hamilton matrix element for a CSF $m$ is given by

[TABLE]

We would like to mention that the $f(b,d)$ and $X(i,j)$ functions are based on and can be expressed as the $D(b,p)$ and $RL$ segment shape function found in Ref. (69). However, we chose the above formulation for conciseness of this manuscript.

B.2 Off-diagonal matrix elements

The off-diagonal matrix element between two CSFs $\ket{m}$ and $\ket{m^{\prime}}$ is given by

[TABLE]

Similar to the Slater-Condon rules Slater (1929); Condon (1930) for matrix element calculation between Slater determinants, we need to identify the involved orbital indices $(i,j,k,l)$ connecting $m^{\prime}$ and $m$ and by comparing the orbital occupation differences between the two CSFs, $\Delta n_{i}=n(d_{i})-n(d_{i}^{\prime})$, already mentioned in Sec. A. There are the following possibilities for $\Delta n_{i}$, which yield a possible non-zero matrix element between $\ket{m^{\prime}}$ and $\ket{m}$:

$\underline{\Delta n_{i}=0}$ for all orbitals, but $\ket{m^{\prime}}$ and $\ket{m}$ differing for some orbitals, implies a full-start $\underline{RL}$ into full-stop $\overline{RL}$ double excitation, type (2c) in Table 6, with only changes in the open-shell orbitals. The matrix element can be expressed as

[TABLE]

Because these full-start into full-stop excitations are symmetric concerning conjugation of the generator indices and the molecular two-body integrals also, it reduces to

[TABLE]

To yield a non-zero matrix element between $m^{\prime}$ and $m$ the indices $i$ and $j$ have to engulf all the differing orbitals, yielding a maximum lower index $I$, and a minimum upper index $J$. Because full-start into full-stop excitations have the possibility to leave $d_{i}$ unchanged, basically all combination $i\leq I$ and $j\geq J$ in the summation have to be considered. There has to be at least two differences between $m$ and $m^{\prime}$, or otherwise it would just be a diagonal matrix element. The singlet coupled $x_{0}$ matrix element branch can be discarded, as a change in $m$ implies $\Delta b=\pm 2$ at least at one orbital. Furthermore the integral, $F(I,J)$, between the region of the first to the last change in $m$ ($I\rightarrow J$) is the same for all matrix elements. The remaining product terms are given by the triplet coupled $x_{1}$ elements for non changing $d_{i}$ value from orbital $i$ to $I$ and $J$ to $j$ given as

[TABLE]

where $RL(d_{i})$ indicates the triplet-coupled $W_{1}(Q_{i};d_{i},d^{\prime}_{i},\Delta b_{i},b_{i})$ matrix elements for a mixed $RL$ generator combination, depending on $d$ and $b$, which can be found in Ref. [(69)]. Since we calculate the matrix elements on-the-fly in the excitation generation step of the FCIQMC method it is useful to formulate the matrix elements in terms of already calculated terms, i.e. $F(I,J)$, to reduce the computational effort of the spin-adapted FCIQMC implementation.

$\underline{\Delta n_{k}=\pm 1:}$ for two spatial orbitals $i,j$. This implies a one-body contribution, as well as two-body contributions with two indices being identical, over which has to be summed. However, there is the additional constraint that the double excitation also has to lead to the same orbital occupancy difference $\Delta n$, which only leaves following terms:

[TABLE]

where the second line, involving weight generators, due to $V_{ijkk}=V_{kkij}$, reduces to

[TABLE]

For $k\neq(i,j)$ both the remaining terms yield (without the two-particle integrals for clarity and $\hat{e}_{ij,kl}=\hat{E}_{ij}\hat{E}_{kl}-\delta_{jk}\hat{E}_{il}$)

[TABLE]

which in total yields

[TABLE]

The third line in Eq. (B.2), due to $V_{ij,kl}=V_{kl,ij}$ and $\hat{e}_{ij,kl}=\hat{e}_{kl,ij}$, reduces to:

[TABLE]

and is a bit more involved to calculate. Depending on the relation of the index $k$ to $(i,j)$, the two-body integral corresponds to certain sequences of generator combinations (assuming $i<j$ for now, which is easily generalized):

$\underline{k<i:}\qquad\underline{LR}\rightarrow\overline{L}R\rightarrow\overline{R}$: type (1j) excitations in Table 6,

without a change in the spin-coupling in the overlap region $(k,i)$ between the two CSFs $m$ and $m^{\prime}$ . The $\Delta b=0$ branch matrix elements of the mixed generator $RL$ contribute multiplicatively $-t^{2}n(d_{k})$ terms, see Ref. [(69)], and the $x_{1}$ matrix element contributions are $x_{1}=\prod_{l=k}^{i-1}RL(l)$, where $RL(l)$ are just the normal mixed generator $x_{1}$ product elements. The product only goes until index $i-1$ to still be able to formulate it in terms of the single excitations $\hat{E}_{ij}$. To formulate it multiplicatively, special factors depending on the step-vector $d(i)$ have to determined, so the semi-stop $x_{1}$ elements of $\overline{R}L(i)/\overline{L}R(i)$ have the same elements as the single starts $\underline{R}/\underline{L}$. The formulation in terms of $\hat{E}_{ij}$ enables us to reuse the already calculated one-body elements in the excitation generation of the FCIQMC method.

The modified values starting $\underline{\Delta}(i)$ and end $\overline{\Delta}(j)$ values for $\underline{RL}\rightarrow R\overline{L}\rightarrow\overline{R}$ and $\underline{RL}\rightarrow\overline{R}L\rightarrow\overline{L}$ type of excitation can be found in Table 14.

The rest of the double excitation overlap matrix elements is the same. So for $k<i$ the two-body matrix elements are given by

[TABLE]

with $\underline{\Delta}(i)$ from Table 14 for the corresponding generator combination.

$\underline{k=i:}\qquad W\underline{R}\rightarrow\overline{R}$: type (0c) excitations in Table 6,

which just reduce to $\bra{m^{\prime}}\hat{E}_{ij}\ket{m}(n(d^{\prime}_{i})-1)$, similarly:

$\underline{k=j:}$

reduces to $\bra{m^{\prime}}\hat{E}_{ij}\ket{m}(n(d_{j})-1)$.

$\underline{k>j:}\qquad\underline{R}\rightarrow R\underline{L}\rightarrow\overline{RL}$: type (1f) excitations in Table 6

As in the $k<i$ case the $x=0$ contribution is a multiplicative $-t^{2}\,n(d_{k})$ factor. The non-vanishing $x_{1}$ overlap matrix elements are again calculated multiplicatively by the use of modified semi-stop segments at $j$ to formulate the matrix element in terms of single excitations for $\underline{R}\rightarrow R\underline{L}\rightarrow\overline{RL}$ and $\underline{L}\rightarrow L\underline{R}\rightarrow\overline{RL}$ generator combinations. The modified terms $\overline{\Delta}(j)$ terms can also be found in Table 14. Here only the step-vector combinations, which lead to the $\Delta b=0$ branch in the overlap region are allowed, since there is not step-vector difference above $j$. The full matrix elements are given by

[TABLE]

with $\overline{\Delta}(j)$ from Table 14 for the specific generator combination.

$\underline{i(j)<k<j(i):}\;\underline{R}\rightarrow\underline{R}\overline{R}\rightarrow\overline{R}\;(\underline{L}\rightarrow\overline{L}\underline{L}\rightarrow\overline{L})$: type (0b) excitations in Table 6

These types of generator combinations correspond to one-body terms actually. At index $k$ the usual product term for the $\hat{E}_{ij}$ matrix element calculation takes on a different than usual value. This modification can be applied multiplicatively, but depends on $d_{k},d^{\prime}_{k},b_{k},\Delta b_{k}$ and the type of generator($i<j$ or $i>j$). The modified values can be found in Table 15.

By defining $r_{k}$ as

[TABLE]

the total matrix element of a single excitation with one-body and two-body contributions can be expressed as

[TABLE]

in terms of the one-body coupling coefficient. The evaluation requires the calculation of the single excitation matrix element $\bra{m^{\prime}}\hat{E}_{ij}\ket{m}$ through Shavitt’s graphical rules Shavitt (1978, 1981), and the summation of terms depending on $d^{\prime}_{k}$ and $d_{k}$ entries of the two CSFs (although in a sequential dependence, since $r_{k}$ depends on $\Delta b_{k}$). We can calculate the $r_{k}$ terms similar to Shavitt’s matrix product terms with $r_{k}$ (88) in the excitation range given in Table 15 during excitation generation. This requires an $\mathcal{O}(N)$ effort in calculation of $r_{k}$, since only occupied orbitals contribute.

$\underline{\Delta n=\pm 1:}$ at two orbitals and additional step-vector differences $\Delta d$ below or above the excitation range, correspond to (1e, 1f, 1i) or (1j) excitations of Table 6, depending on the ordering of the remaining indices. These are $d=1,d^{\prime}=2$, and vice versa, step-vector differences outside the range $(i,j)$, corresponding to double excitations with mixed generator full-starts $\underline{RL}$ or full-stops $\overline{RL}$. Similar to the $\Sigma\lvert\Delta n\rvert=0$ case all possible excitations connecting the two CSFs have to engulf the first step-vector change, $I$, if it occurs before $\min(i,j)$ or the last step-vector change, $J$, if it is after $\max(i,j)$. However, all mixed full-starts before $I$ or after $J$ have to be considered too, since there is the possibility of a $\underline{RL}(d=1,d^{\prime}=1)$ or $\underline{RL}(d=2,d^{\prime}=2)$ start with non-zero $x_{1}$ matrix element Shavitt (1981). So the matrix element is given by:

[TABLE]

with $RL(d_{l})$ again being the $x_{1}$ matrix elements, $F(I,i)/F(j,J)$ being the always involved $x_{1}$ matrix elements engulfing all step-vector changes in the overlap region and $\underline{\Delta}(k)/\overline{\Delta}(k)$ being the modifying terms to express it in terms of single excitation matrix elements $\hat{E}_{ij}$, see Table 14.

$\underline{\Delta n_{i}=\pm 2}$ at two spatial orbital $i$ and $j$ implies a full-start into full-stop double excitation with two alike generators ($\underline{RR}\rightarrow\overline{RR}$ or $\underline{LL}\rightarrow\overline{LL}$, corresponding to type (2a) and (2b) in Table 6). This completely specifies the indices and the full matrix element is just given by

[TABLE]

and calculated with Shavitt’s graphical rules Shavitt (1981). The order of the orbitals, where electrons are removed and added, determines the type of generators.

$\Delta n\neq 0$ at three different orbitals with one $\Delta n=\pm 2$ and two $\Delta n=\mp 1$, corresponds to type (1c, 1d, 1g) or (1h) excitations of Table 6. This determines all four indices, with two indices being identical, with $\Delta n_{k}=\pm 2$ and $\Delta n_{i}=\Delta n_{j}=\mp 1$. Leaving the matrix element to be:

[TABLE]

$\underline{\Delta n_{i}\neq 0}$ at four different orbitals with two times $\Delta n=1$ and two times $\Delta n=-1$ values, corresponds to type (3a, 3b, 3c, 3d, 3e) or (3e) excitations of Table 6. This also completely determines all four indices of the excitation, but there are four different combinations of these indices which can lead to the same state, where two of them, however, are equivalent. The relation and ordering of these indices determines the type and combinations of generators, with the total matrix element given by

[TABLE]

where for orbitals $l$ and $k$, $\Delta n_{l}=\Delta n_{k}=+1$ and for $i$ and $j$, $\Delta n_{i}=\Delta n_{j}=-1$. The relative positions of the $\Delta n=+1$ and $\Delta n=-1$ orbitals determines the generator combinations and type of excitations involved.

For alike generator combinations, e.g. $\underline{R}\rightarrow\underline{R}R\rightarrow\overline{R}R\rightarrow\overline{R}$, we have to take into account the sign flip due to an exchange of operator indices. Since $\hat{e}_{ik,jl}$ and $\hat{e}_{il,jk}$ ($i<j<k<l$), both contribute to the same excitation. As already pointed out by Paldus Paldus and Boyle (1980); Boyle and Paldus (1980), the Coulomb and exchange type contributions can be expressed in terms of the same $x=0$ and $x=1$ matrix element contributions with

[TABLE]

where $w_{0}\neq 0$ only if $\Delta b_{k}=0,\forall k\in S_{1}$. By sticking to the convention to use the standard order of operators, as indicated in Table 6, the contribution of an exchange of orbital indices in the generator can be expressed as

[TABLE]

with a type (3a) excitation from Table 6 as an example. The total matrix element is then given by

[TABLE]

The case of alternating orbital occupancy differences and $\Delta n=\pm 1\rightarrow\Delta n=\mp 1$ involve no sign change in mixed generator semi-start and semi-stops for the $x_{1}$ matrix element. These type (3c-3f) excitations of Table 6 also contain non-overlap generator combinations, indicated by the subscript 0. Because the non-overlap double excitations are contained as the $\Delta b_{k}=0,\forall k\in S_{1}$ special case of these mixed generator excitation, we do not treat them explicitly, but stick to the convention to always use the mixed generator combination in the excitation generation. For an excitation, which left the $\Delta b_{k}=0$ path at some point in the overlap region, only the $x=1$ matrix element contributes. If $\Delta b_{k}=0,\forall k\in S_{1}$, the Coulomb type contribution can be obtained by the $x=0$ term of the exchange matrix element Paldus and Boyle (1980)

[TABLE]

otherwise $w_{0}=0$. The total matrix element is then given by

[TABLE]

$\underline{\Delta n\neq 0}$ at more than four different spatial orbitals or $\underline{\Delta b\neq 0}$ outside of excitation range for $\Sigma\lvert\Delta n\rvert=4$, yields a zero matrix element, as such excitations cannot be obtained by a single application of the Hamiltonian.

For brevity of this manuscript the extensively used remaining two-body segment value tables are not listed here; the interested reader is referred to References [(43; 69)] for a detailed explanation and listing of the GUGA matrix elements by I. Shavitt.

Appendix C Weighted Orbital Choice with GUGA Restrictions

A note on the weighting of the integral contribution to the orbital picking process: As one can see in Eq. (97) and (100) it is not as easy as in a SD based implementation to weight an integral contribution of orbitals by the exact matrix element, $V_{ikjl}$ for spin-opposite and $V_{ikjl}-V_{iklj}$ for spin parallel excitations, as it is done in the current FCIQMC implementation. Since the relative sign of the $w_{0}$ and $w_{1}$ contribution (95, 96) depends on the chosen excitation $\ket{m^{\prime}}$ and can not be easily predetermined and there is no notion of a $m_{s}$ quantum number in a spin-adapted calculation. Our choice was to weight the integral contribution by the magnitude of the integrals $\lvert V_{ikjl}\rvert+\lvert V_{iklj}\rvert$ to capture the strongest couplings at least. This leads to some inefficiencies in the excitation generation of the CSF based implementation.

C.1 Restrictions on the Orbital Choice for Single Excitations

To ensure at least one possible non-zero excitation, $\hat{E}_{ij}\ket{m}$, we have to place some additional restriction on the choice of orbital $(i,j)$ compared to a SD-based implementation. The idea is to first pick an electron in an occupied spatial orbital $j$ at random with $p(j)=1/N$. Depending on the step-value $d_{j}$ certain restriction on the to-be-picked orbital $i$ are placed. A general restriction is that $i$ must not be doubly occupied $d_{i}\neq 3$.

If $d_{j}=3$, since both $\Delta b$ branches can end at $d_{j}$ for raising generators R $(i<j)$ and also both branches can start for a lowering generator L $(i>j)$ there are no additional restrictions on the orbital $i$, except $d_{i}\neq 3$.

For $d_{j}=1$ there is only a $\Delta b_{j}=+1$ start $\underline{L}$ and a $\Delta b_{j-1}=-1$ end $\overline{R}$ possible. So there is the restriction, that $d_{i}$ must be [math], which allows both $\Delta b$ branches to start or to end, or $d_{i}=2$, which would lead to the correct $\Delta b_{i}=-1$ start for $\underline{R}$ or would allow the $\Delta b_{i-1}=+1$ end for $\overline{L}$. A $d_{i}=1$ value is only allowed, if there is a valid switch possibility $d_{k}=2$ in the range $(i,j)$. For a chosen $d_{j}=2$ electron orbital the restrictions are similar with $d_{i}=\{0,1\}$ being valid in general, and $d_{i}=2$ only if a switch possibility $d_{k}=1$ for $k\in(i,j)$. The actual restriction is implemented by finding the adjacent opposite spin-coupled orbitals $i_{lower}$ and $i_{upper}$ for a $d_{j}=\{1,2\}$ and only allowing $d_{i}=d_{j}$ to be picked if $i<i_{lower}$ or $i>i_{upper}$. A flow-chart of this decision-making process is given in Fig. 20.

If we want to make use of an Abelian point group symmetry, e.g. $D_{2h}$, it is much easier than suggested in the literature Brooks et al. (1980), to just restrict the choice of orbital $i$ to the symmetry allowed orbitals $n_{j}$ corresponding to the picked electron orbital $j$, since for single excitations, in order for $t_{ij}$ to be nonzero, the product of the irreducible representations for orbitals $i$ and $j$ must be totally symmetric.

C.2 Restriction on the Orbital Choice for Double Excitations

The excitation generation for doubles is a bit more involved, but due to the product structure of the matrix elements (72). We follow the same approach as for single excitations to pick the four orbitals of $\hat{e}_{ij,kl}$ in such a manner to have the probability $p(ijkl)$ be related to the integral contribution of the Hamiltonian matrix element and at the same time ensure that at least one valid excitation can be reached.

We first pick the ordered electron pair $(j<l)$ at random with a probability $p(jl)=1/N_{pairs}$, where $N_{pairs}$ is the number of electron pairs in the simulation. The first orbital $(i)$ to excite to is then picked out of all, non-doubly occupied $d_{i}\neq 3$ orbitals, weighted with the Cauchy-Schwarz inequality based approximation Smart, Booth, and Alavi ; Holmes, Changlani, and Umrigar (2016) of the integral contribution $V_{ijkl}$.

The major change, compared to a SD based implementation, now comes only in the choice of the second orbital $(k)$ to excite to. Here we place the restrictions depending on the possible used spatial symmetry and the additional restriction, due to the UGA to obtain non-zero excitations. Additionally, we restrict the picking of orbital $(k)$ in such a way that we do not pick *quasi-*single excitations, which are already taken account of in the single excitation matrix elements. The overall restriction $d_{k}\neq 3$ remains of course.

To formulate the conditions for a valid orbital index choice, $(i,j,k,l)$, we have to look at the properties of the non-zero two-body segment shapes. The semi-start segments, $\underline{L}L,\underline{R}R$, behave similar to single segment shapes concerning the restrictions on the in-coming and out-going $\Delta b$ values of an excitation and are listed in Table 16. And similar to the end of a single excitation there are certain restriction for non-zero two-body elements at the end of the overlap range, depending on the type of the two alike generators, see Table 17. As one can see in these tables these segments behave like a single-excitation starts for an in-going $\Delta b_{k-1}=0$ branch and like a single excitation end segment for the approaching $\Delta b_{k-1}=\pm 2$ branches.

For two-body generators with identical starting indices, $\hat{e}_{ij,ik}$, only the $\Delta b_{k}=0$ branch contributes, due to a zero $x=1$ matrix element in the overlap region Shavitt (1981). In addition only a $\Delta b_{k}=0$ branch leads to a non-zero matrix element with two coinciding upper indices, $\hat{e}_{ij,kj}$. This means that for the $\underline{R}R$ or $\underline{L}L$ segments only the $\Delta b_{k}=0$ branch can be chosen. So these type of double excitations can be treated very similar to single excitations, since the $x=0$ contribution is very easy to compute, ($-1^{n_{2}}$, with $n_{1}$ being the number of singly occupied orbitals in the overlap region) and are only non-zero if $d_{k}^{\prime}=d_{k}$ in the overlap range. So no switch decisions have to be made in the excitation generation.

The case of mixed generators $R+L$ is a bit more involved. A simultaneous start $e_{ij,ki}$ acts similar to an usual double intermediate segment value, except that a $d_{i}=0$ value leads to a zero matrix element, see Tables 18 and Ref. [(69)]. And similar to intermediate segment values of alike generators in the overlap region, the $x=0$ matrix elements are zero for the $\Delta b_{k}=\pm 2$ branches.

The $x=1$ contribution of the mixed two-body segment values, see Ref. [(69)], is zero for $d_{i}=3$, but not for the $\Delta b_{k}=0$ branches of $d_{i}=\{1,2\}$. This leads to a major complication in the implementation of CSFs in the FCIQMC algorithm through the GUGA approach. These contributions with no change in step-value with a non-zero matrix element correspond to an exchange type contribution to double excitations. Since in the FCIQMC excitation generation it is necessary to uniquely assign a definite probability $p(m^{\prime}|m)$, different starting orbitals $i^{\prime}<I$, with $I$ indicating the first step-value change $\Delta d_{I}$, can contribute to an excitation with a $\underline{RL}$ start. The matrix element influence was mentioned above, but also the probabilities, $p(i^{\prime})$, for all possible other starting orbitals $i^{\prime}<I$ have to be accounted for. Similarly for an $\overline{RL}$ end, all other possible $j^{\prime}>J$, with $J$ indicating the last step-value change $\Delta d_{J}$, ending orbitals have to taken into account. And for a pure exchange type excitation (type 2c in Table 6) $\underline{RL}\rightarrow\overline{RL}$ all combinations $(i^{\prime}<I,j^{\prime}>J)$ of possibly contributing orbitals have to be considered in the matrix element and generation probability computation.

Otherwise a $R\underline{L}$ segment behaves similar to an $L\underline{L}$ and $\underline{R}L$ to an $\underline{R}R$ in terms of $d_{k}^{\prime},b_{k}$ and $\Delta b_{k-1}$ restrictions, except the number of electrons in $N_{k}^{\prime}=N_{k}$. Also the intermediate $RL$ segments behave as $LL$ and $RR$ and the $R\overline{L}$ is equivalent to $L\overline{L}$ and $L\overline{R}$ to $R\overline{R}$ respectively. Except the electron number difference becomes the corresponding value $N_{k}^{\prime}=N_{k}\pm 1$ of the ongoing excitation ($R$ in the case of $R\overline{L}$ and $L$ for $L\overline{R}$). It should also be noted, that a $\Delta b_{k}=0$ branch can end at any $d_{j}\neq 0$ value, whereas $\Delta b_{k}=-2$ is restricted to $d_{j}=1$ and $\Delta b_{k}=+2$ to $d_{j}=2$, to be able to align the $S_{k}$ value of $\ket{m^{\prime}}$ and $\ket{m}$, so they coincide outside of the range of the generator $\hat{e}_{ij,jl}$. Since $\Delta b_{k}=0$ already indicates $\ket{m^{\prime}}=\ket{m}$ in the overlap range, this issue is no problem in a direct CI calculation with CSFs, but is burdensome to implement in FCIQMC, since we want to be able to get one out of all possible excitations for a given CSF $\ket{m}$ and assign a unique generation probability to it. So we also have to take into account all other possible index combinations, which would be able to lead to this excitation and sum their matrix elements of course, but also recompute the probability.

Due to the uniqueness for most type of excitations, this is no problem, except in the case of these exchange type excitations with coinciding indices and a $RL$ generator combination. Unfortunately we have not yet found a more elegant way to treat these cases, except implement it in the most efficient way. With a heavy re-usage of terms to avoid an $\mathcal{O}(N^{2})$ or even $\mathcal{O}(n^{2})$ computational cost of these excitations.

C.3 Orbital Picking and Excitation Identification

In the following, the work flow of picking a valid index combination $(i,j,k,l)$ for a non-zero double excitation of a CSF $\ket{m}$ in the FCIQMC method is presented. A flow-chart of the decision-making process is shown in Fig. 21.

Both electron indices $(j<l)$ and the first hole index $(i)$ are picked with uniform or a weighted probability. If both picked electrons are in the same spatial orbital $\underline{j=l\rightarrow d_{j}=3}$, similar restrictions as for single excitations apply, for the remaining orbitals $i$ and $k$. Of course both orbitals $i$ and $k$ have to be non-doubly occupied. This applies in general, independent of the relation of electron orbitals $j$ and $l$ and their step-value $d_{j},d_{l}$. And $i=k$ is only possible if $d_{k}=0$, since both electrons will be excited to the same spatial orbital. Since both electrons get removed from the same orbital the only possible excitation types are (1b,1c,1d,2a,2b) of Table 6. If $i=k$ the type of excitation is (2a) if $i>j$, or (2b) if $i<j$, requiring $d_{i}=0$. For $i\neq k$ the same restrictions as for single excitations apply, that $d_{i}=d_{k}=1$ is only possible if a switch possibility $d_{m}=2$ in the range $(i,k)$ and vice versa for $d_{i}=d_{k}=2$. The type of excitation only depends on the order of the involved indices. As already mentioned, all these excitations require $\Delta b_{m}=0$ in the overlap range. With the trivial case of type (1b) excitation with a single orbital overlap range. Which make the calculation of the excitation very similar to single excitations.

If the electron indices are not equal $j\neq l$, the picking of the remaining orbitals $k$ and $i$ depends more strongly on the step-values of the already chosen orbitals $(j,l)$.

If $\underline{d_{j}=d_{l}=3}$ there is no additional restriction on orbital $k$, since in except of $b$ value restrictions on starts of a doubly occupied orbital ($\Delta b=+1$ branch forbidden due to $b=0$, e.g.) all restrictions mentioned in the previous section can be accounted for, due to the flexibility of the $d_{j}=d_{l}=3$ step-values.

If $i=k$, depending on the order of the indices this leads to excitations of type (1a,1d) or (1g), since $i\neq j,l$ due to $d_{j}=d_{l}=3$. These excitations again can be easily treated, due to the single overlap region (1a) or a necessary $\Delta b_{m}=0$ in the overlap region.

If $i\neq k$ all 4 indices are different, leading to a type (3*) excitation depending on the order of the indices. Where again, the exchange type excitation is chosen by definition and not a possibly non-overlap double excitation (3c0,3d0,3e0,3f0).

If $\underline{d_{j}=3,d_{l}=\{1,2\}}$: Depending on if the already picked first orbital to excite to $i=l$, orbital $k$ must be restricted to $k<l$. This is because $k>l$ would lead to an exchange contribution to a single excitation, which is already taken into account for in the singles matrix elements. This leads to excitations of type (1e) or (1f) depending on the order of $k$ and $j$. If $i>l$ orbital $k\neq l$ since this again would lead to an already accounted exchange contribution to a single excitation. If $i=k$ it is an (1d) excitation otherwise it is one of the type (3*) depending on the relation of the indices. If $i<l$ there is no restriction on the indices for $k$ and this can lead to a variety of excitations.

If $\underline{d_{j}=\{1,2\},d_{l}=3:}$ There are similar restrictions considering the already picked orbitals. If $i=j$, $k$ must be picked $k>j$ to avoid choosing already accounted for exchange contributions to single excitations. And if $i<j$ $k$ must not be $j$, otherwise there are no additional restrictions.

If both $\underline{d_{j}=\{1,2\}}$ and $\underline{d_{l}=\{1,2\}}$ are singly occupied the most stringent restrictions apply. If $i=l$ there should be a possible switch between the already picked $j$ and $l$ if both have the same step-value $d_{j}=d_{l}$, since otherwise it would not be possible to fulfill the $\Delta b$ criteria at the end of an excitation to lead to a non-zero excitation. Additionally, orbital $k$ has to be lower than $i$, otherwise it is again an exchange contribution to a single excitation. Depending on the order of the orbitals, this leads to type (1e,1f,2c) excitation, which, already mentioned, needs additional re-computation of matrix element and generation probability contribution.

If $i=j$, there also must be a switch possibility for $d_{j}=d_{l}$ between $j$ and $l$ and $k>j$. This leads to type (1i,1j,2c) excitations. With the necessity of recalculation of matrix element and generation probability contributions.

If $i>j$ orbital $k$ must not be $l$ to avoid an exchange contribution to singles and similarly if $i<j$, $k$ must not be $j$. And $k$ can only coincide with $j$ for $i>j$, if there is a switch possibility between $j$ and $l$, if $d_{j}=d_{l}$. And similar for $i<j$, $k=l$ is only possible if $d_{j}\neq d_{l}$ or there is a switch possibility between $j$ and $l$. Otherwise no restrictions are place on the picking of orbital $k$ and the type of excitation depends on the order of the indices and can lead to all sort of excitation types in Table 6.

In the whole picking process, since we allow the empty orbitals to be picked in any order, in addition to $p(i|jkl)$ we also have to recompute $p(k|ijl)$ of having picked the orbitals in the opposite order, since they lead to the same possible excitation. This increases the generation probability by 2 in general, but introduces the effort to recompute. We could, similar to the electron orbitals $j$ and $l$, decide to pick only orbitals $i<l$, which would also make the identification of the excitation type easier.

To utilize Abelian point group symmetry for double excitations we restrict the last picked empty orbital $k$ to ensure that the product of the four irreducible representations of $(i,j,k,l)$ are totally symmetric and consequently $V_{ijkl}$ is not zero due to symmetry.