Singlet-Triplet Energy Gaps of Organic Biradicals and Polyacenes with Auxiliary-Field Quantum Monte Carlo

TL;DR

This paper demonstrates that phaseless auxiliary-field quantum Monte Carlo (ph-AFQMC) can accurately predict singlet-triplet energy gaps in biradical molecules and large polyacenes, enabling reliable modeling of photochemical processes.

Contribution

The study introduces a scalable ph-AFQMC methodology with optimized trial wavefunctions for accurate singlet-triplet gap predictions in complex biradical systems and large polyacenes.

Findings

Achieves ~1 kcal/mol accuracy in predicting ST gaps for small molecules.

Validates the method's effectiveness on large polyacenes like pentacene.

Provides a protocol for selecting trial wavefunctions based on spin-contamination.

Abstract

The energy gap between the lowest-lying singlet and triplet states is an important quantity in chemical photocatalysis, with relevant applications ranging from triplet fusion in optical upconversion to the design of organic light-emitting devices. The ab initio prediction of singlet-triplet (ST) gaps is challenging due to the potentially biradical nature of the involved states, combined with the potentially large size of relevant molecules. In this work, we show that phaseless auxiliary-field quantum Monte Carlo (ph-AFQMC) can accurately predict ST gaps for chemical systems with singlet states of highly biradical nature, including a set of 13 small molecules and the ortho-, meta-, and para-isomers of benzyne. With respect to gas-phase experiments, ph-AFQMC using CASSCF trial wavefunctions achieves a mean averaged error of ~1 kcal/mol. Furthermore, we find that in the context of a spin-projection technique, ph-AFQMC using unrestricted single-determinant trial wavefunctions, which can be readily obtained for even very large systems, produces equivalently high accuracy. We proceed to show that this scalable methodology is capable of yielding accurate ST gaps for all linear polyacenes for which experimental measurements exist, i.e. naphthalene, anthracene, tetracene, and pentacene. Our results suggest a protocol for selecting either unrestricted Hartree-Fock or Kohn-Sham orbitals for the single-determinant trial wavefunction, based on the extent of spin-contamination. These findings pave the way for future investigations of specific photochemical processes involving large molecules with substantial biradical character.

Figures

Tables

Table 1. Table 1: ⟨ S 2 ⟩ delimited-⟨⟩ superscript 𝑆 2 \langle S^{2}\rangle of the singlet and triplet unrestricted solutions for the 13 small molecule biradical set.

	UHF singlet	UHF triplet	UB3LYP singlet	UB3LYP triplet
CH2	0	2.02	0.57	2.01
CF2	0.14	2.01	0	2.00
NH${}_{}{}^{+}{}_2{}^{}$	0.82	2.02	0.75	2.01
SiH2	0.43	2.01	0	2.00
PH${}_{}{}^{+}{}_2{}^{}$	0.44	2.01	0	2.00
H2CC	0.16	2.35	0	2.03
NH	1.01	2.02	1.00	2.01
NF	1.01	2.02	1.00	2.01
O2	1.02	2.05	1.01	2.01
OH+	1.01	2.01	1.00	2.01
O	1.01	2.01	1.00	2.00
C	1.02	2.01	1.01	2.00
Si	1.05	2.02	1.01	2.00

Table 2. Table 2: Mean signed, absolute, and maximum errors [kcal/mol] of the theoretical methods shown in Fig. 1 for the 13 small molecule biradical set. Sorted by MAE value. Parenthesis indicate statistical errors.

	MSE	MAE	MaxE
AFQMC/CAS	0.5(4)	0.9(4)	1.7(3)
AFQMC/U	0.5(7)	1.2(7)	2.1(8)
W2X	2.7	3.1	6.9
WA-KS/BLYP	-6.8	7.3	17.4
UKS/BLYP	-12.5	12.9	34

Table 3. Table 3: Mean signed, absolute, and maximum errors [kcal/mol] of ph-AFQMC results and other methods for the 8 molecule biradical subset shown in Fig. 2 . Sorted by MAE value. Parenthesis indicate statistical errors.

	MSE	MAE	MaxE
AFQMC/U$†$	0.83(7)	1.0(7)	2.1(7)
AFQMC/CAS$†$	0.5(4)	1.1(4)	1.7(3)
SF-LDAa	-0.7	1.7	4.0
SF-CIS(D)b	1.8	1.9	5.1
W2Xc	3.0	3.7	6.9
(V)FS-PBEd	4.2	4.3	6.7
B3LYP/pp-RPAe	-3.5	4.8	7.7
WA-KS/BLYPc	-6.6	6.6	17.5
HF/pp-RPAe	-7.0	7.1	17.3
UKS/BLYPc	-12.5	12.5	34.0

Table 4. Table 4: ⟨ S 2 ⟩ delimited-⟨⟩ superscript 𝑆 2 \langle S^{2}\rangle of the singlet and triplet unrestricted solutions for ortho- meta- and para- benzyne molecules.

	UHF singlet	UHF triplet	UB3LYP singlet	UB3LYP triplet
ortho	1.26	2.31	0.00	2.01
meta	0.97	2.68	0.10	2.02
para	1.68	2.31	0.92	2.01

Table 5. Table 5: ST gaps [kcal/mol] for the ortho- meta- and para- benzyne isomers. Parenthesis indicate statistical errors.

	ortho	meta	para
expt*	37.5 $\pm$ 0.3	21.0 $\pm$ 0.3	3.8 $\pm$ 0.3
ZPE**	-0.6	1.0	0.5
ZPE-corr’d expt	38.1	20.0	3.3
AFQMC/CAS$†$	37.4(6)	20.7(8)	4.5(5)
AFQMC/U$†$	37.6(7)	18.9(9)	2.2(9)
UB3LYPa	29.4	14.2	2.4
HF/pp-RPAb	45.6	35.5	4.0
B3LYP/pp-RPAb	37.4	22.1	0.6
SF-CIS(D)c	35.7	19.4	2.1
SF-B3LYPd	46.9	26.1	6.9
SF-CCSD(T)e	37.3	20.6	4.0
SF-oo-CCDc	37.6	19.3	3.9

Table 6. Table 6: ⟨ S 2 ⟩ delimited-⟨⟩ superscript 𝑆 2 \langle S^{2}\rangle of the singlet and triplet unrestricted solutions for polyacenes n = 2 − 5 𝑛 2 5 n=2-5 .

$n$	UHF singlet	UHF triplet	UB3LYP singlet	UB3LYP triplet
2	1.10	2.30	0	2.02
3	1.78	2.68	0	2.02
4	2.43	2.91	0	2.03
5	3.06	3.44	0	2.03

Table 7. Table 7: ST gaps [kcal/mol] for the polyacenes n 𝑛 n =2-5. Square brackets indicate alternate experiments, and parenthesis indicate statistical errors.

	naphthalene	anthracene	tetracene	pentacene
expt*	[60.9] 61.0	[42.6] 43.1	29.4	19.8$\pm$0.7
ZPE*	-3.4	-2.3	-1.8	-1.5
ZPE-corr’d expt	64.4	45.4	31.2	21.3
AFQMC/Ua	68.0(1.2)	46.2(1.2)	34.0(1.6)	25.2(1.6)
UB3LYPb	62.6	41.8	27.7	17.9
CCSD(T)/FPAc	65.8	48.2	33.5	25.3
B3LYP/pp-RPAd	66.2	45.7	32.1	22.6
GAS-pDFT (FP-1)e	70.6	45.5	33.6	25.4
GAS-pDFT (WFP-3)e	64.7	43.1	28.8	20.5
ACI-DSRG-MRPT2f	62.2	43.2	28.3	18.0
DMRG-pDFTg	67.1	46.1	31.6	22.6

Full text

\alsoaffiliation

Department of Physics, College of William and Mary, Williamsburg, VA 23187

1 Introduction

The energy gap separating the lowest-lying singlet and triplet states of a molecule is an important property relevant to many chemical processes. For example, light absorption by chlorophyll in Photosystem II can produce triplet states which in turn react with triplet oxygen to produce short-lived and highly reactive singlet oxygen.1 Additionally, the relative energetics of first-excited singlet and triplet states in dopants utilized in organic light-emitting diodes governs the efficiency of such devices, and is a useful parameter for the design of light-emitting electronics2. In the field of photocatalysis the singlet-triplet (ST) gap is directly relevant to a variety of redox reactions.3 In addition the ST gap is a quantity of crucial importance for determining the energetic feasibility of the optical processes known as singlet fission4 and upconversion5. In the former process, a single photon produces a high energy singlet state which eventually splits into two triplet excitons; in the latter, two triplet excitons annihilate or fuse to form a high energy emissive singlet. In certain cases the ST gap can be probed under cryogenic conditions via phosphorescence measurements, however for a large number of relevant molecules direct measurement is not possible.

The electronic structures relevant to these types of applications can be complicated by the presence of biradical character in one or more of the involved spin-states. Biradicals are molecules in which two valence electrons can occupy two degenerate but spatially distinct molecular orbitals6 (the species are referred to as biradicaloids if these two orbitals are nearly-degenerate, but we do not make this distinction here). The many possible electronic configurations give rise to their capacity to exhibit remarkably specific chemical reactivity.7, 8, 9, 10, 11. Singlet states can be characterized as either closed-shell or open-shell. In the former, one of the two valence states is doubly occupied, which is typically the case, e.g., in carbenes. These species can simultaneously display Lewis base and Lewis acid character, and thereby undergo concerted addition reactions, e.g. with alkenes, to produce stereospecific products. Open-shell singlet states12 are characterized by single-occupancy of each of the two valence states, as in all triplet states due to Pauli exclusion. Such open-shell molecules yield products of mixed stereochemistry.

The electronic structure community has witnessed the development of promising theoretical methods for computing the ST gap in biradical molecules. This is a challenging problem for traditional single-reference ab initio methods, e.g. Hartree-Fock (HF) and Density Functional Theory (DFT), as a minimal quantum-mechanical description of the wavefunctions corresponding to all singlet states and one triplet state ($M_{s}=0$) in the two electron two orbital model is necessarily a superposition of two electronic configurations,13, 14 thus requiring more than one Slater determinant. Even for simple chemical species such as NH and O2, which have triplet ground states, the singlet spin-configuration is notoriously difficult to describe. Accurate predictions are further complicated by the requirement that static and dynamic electron correlation be well balanced, with spin-states of both multiplicities treated on equal theoretical footing. Biradical systems have been studied with DFT15, 16 and fractional spin variants,17, 18 Generalized Valence Bond theory,19, 20, 21, 22 spin-flip (SF) methods,13, 23, 24, 25 Complete Active Space Self-Consistent Field (CASSCF) with 2nd order perturbation theory,26 multiconfiguration pair-density functional theory (pDFT),27, 28, 29, 30, 31 Coupled Cluster (CC) methods,32, 33, 34 doubly electron-attached equation-of-motion CC theory, 35, 36 spin-extended Configuration Interaction (CI) with singles and doubles, 37 incremental Full CI,38 Difference Dedicated CI,39, 40, 41 the particle-particle random phase approximation (pp-RPA)14, 42 the Density Matrix Renormalization Group (DMRG),43 Yamaguchi spin projection44 and its recent combination with orbital-optimized MP245.

In this work we use phaseless auxiliary-field Quantum Monte Carlo (ph-AFQMC)46, 47 to compute ST gaps. While imaginary-time projection is most frequently used to yield ground-state properties of a system, the formalism has also been used to accurately compute low-lying excited states of materials48, molecular diatomics49 and dipole-bound species.50 These calculations rely on the fact that eigenfunctions of the Hamiltonian are orthogonal. In practice, when the exact eigenfunctions are unknown, ph-AFQMC calculations use a so-called “trial wavefunction” to project out orthogonal components that may be sampled along the imaginary-time random walks. For example, for the molecules with triplet ground-states relevant to this work, the energy of the lowest-lying excited singlet state can be sampled using a trial wavefunction with $\langle S^{2}\rangle=0$, due to its near-orthogonality to the true triplet ground-state. In the limit of an exact trial wavefunction, this symmetry-constraining approach would be exact.

In what follows we will show that an unrestricted single-determinant trial wavefunction is capable of accurately describing multi-reference biradical species. This is a significant result because the computational cost of CI and CASSCF calculations, despite many recent advances,51, 52, 53, 54 scales exponentially with system size and thus such methods are infeasible as trial wavefunctions for ph-AFQMC. As an example, the number of $\pi$ electrons in the polyacene series is 4$n$+2 ($n$=1,2,… for benzene, naphthalene,…), which typically must be included in the active space for accurate MCSCF-based predictions. Typical CI solvers can handle up to $\sim 16$ active electrons and orbitals, which would be insufficient for $n>3$.

We view the ability of an electronic structure theory to accurately describe biradicals as a prerequisite for future studies of large, typically conjugated systems that catalyze photochemical processes such as upconversion. After showing that a spin-projected approach to ph-AFQMC55 with unrestricted single-determinant trial wavefunctions produces accurate results for a set of strongly biradical small molecules and three benzyne isomers, we then illustrate the scalability of our approach by taking a first step toward relevant and relatively large photocatalytic molecules, namely the polyacene series including naphthalene, anthracene, tetracene, and pentacene. These molecules are well-studied both experimentally (with measurements available in the literature for $n$=1-5, though not beyond) and theoretically, as they have myriad applications in organic electronics (see references in the first paragraph of Ref. 42). While biradical and polyradical character is predicted to be responsible for the instability of hexacene and longer acenes,43, 56 we focus on the $n$=2-5 molecules since they are representative of the majority of molecules in our target class of photocatalysts.3 In fact, $n$=3,4 are known to perform upconversion,5 and $n$=3-5 are known singlet-fission catalysts.4

This work is organized as follows. Details of our computational approach are described in Sec. II. In Sec. III we compute ST gaps for a set of 13 small organic molecules which have singlet states of highly biradical nature, and compare with experimental measurements. Next we examine ortho- meta- and para- isomers of benzyne, and show that very high accuracy can be obtained with both CASSCF and single-determinant trial wavefunctions using a basis set of moderate size. Having shown that a single-determinant trial wavefunction combined with a spin projection technique is capable of accurately describing multi-reference biradical species, we proceed to compute ST gaps with ph-AFQMC for the increasingly large (but not necessarily biradical) systems naphthalene, anthracene, tetracene, and pentacene, and compare our results with state-of-the-art electronic structure theories and experimental measurements.

2 Computational Details

The ph-AFQMC methodology is reviewed in Refs. 57 and 58. All calculations utilize an imaginary time step of $\Delta\tau=0.005Ha^{-1}$. Walker orthonormalization, population control, and local energy measurements are carried out every 2, 20, and 20 steps, respectively. We utilize a modified Cholesky decomposition of the electron repulsion integrals,59 with cutoffs of $10^{-5}$ for the small molecule biradicals, and $10^{-4}$ for polyacenes, $n$=2-4. For pentacene ($n$=5) we use density fitting with the Weigend Coulomb-fitting basis set60 to reduce the memory requirements of the calculation, while preserving high accuracy in energy differences61.

All calculations utilize single-precision floating point arithmetic. For ionization and bond-dissociation energies of transition metal atoms and diatomics, respectively, this yielded very high accuracy while reducing the computational cost compared to double-precision calculations61, 62. For pentacene, the largest molecule considered in this work, we verified with separate calculations in the STO-3G basis that single- and double-precision calculations gave statistically indistinguishable results.

For the small molecule biradicals, we use the aug-cc-pV$x$Z basis sets,63 with $x$=T,Q. Unrestricted HF (UHF), restricted HF (RHF) and its open-shell variant (ROHF), unrestricted Kohn-Sham DFT with the B3LYP functional (UB3LYP), and CASSCF trial wavefunctions are obtained using PySCF64. ST gaps are extrapolated to the complete basis set (CBS) limit using exponential and $1/x^{3}$ forms for the mean-field and correlation energies, respectively, as detailed in, e.g., Ref. 61.

For the benzyne isomers and polyacenes we report ST gaps in the cc-pVTZ basis,65 primarily for computational expedience, since pentacene has nearly 900 basis functions. For hydrocarbon systems, full CBS extrapolation typically alters the triple-zeta results by $\sim$1 kcal/mol or less. This empirical finding is consistent with previous studies of polyacenes using other wavefunction methods which show little basis set dependence.32, 66, 42

Our ph-AFQMC calculations utilize the spin-projection technique detailed in Ref.55. The walkers are initialized with RHF for singlets and ROHF for triplets, which have $\langle S^{2}\rangle$ of exactly 0 and 2, respectively. With an appropriate form of the Hubbard-Stratonovich transformation this can ensure that the single-particle imaginary-time propagator preserves spin-symmetry despite the use of a (possibly) spin-contaminated unrestricted trial wavefunction.

In what follows, we propose a simple empirical protocol, “AFQMC/U”, to guide the selection of an optimal unrestricted orbital set, among UHF and UB3LYP orbitals, to be used as the trial wavefunction in AFQMC. UHF orbitals are used by default, unless $\langle S^{2}\rangle_{singlet}^{UHF}$ or $\langle S^{2}\rangle_{triplet}^{UHF}$ deviate by more than 10$\%$67 from their spin-pure values, i.e. 0 and 2, respectively. For biradical singlet states, UB3LYP orbitals are utilized only when $\langle S^{2}\rangle_{singlet}^{UHF}>1.1$, since wavefunctions constructed from unrestricted Kohn-Sham orbitals are known to exhibit less spin contamination compared to UHF solutions.68 Biradical triplet states are not expected to exhibit significant spin-contamination as they are typically well-described by a single determinant (the two-determinant $M_{s}$ = 0 triplet state is not encountered due to the constraint on $\langle S_{z}\rangle$ in the quantum-chemical methods relevant here).

3 Results and Discussion

3.1 Small Molecule Biradical Set

ST gaps are computed for the set of 13 small molecules with singlet states that exhibit substantial biradical character, recently examined in Ref. 27. Molecular geometries for the biradical species are taken from the Supporting Information of Ref. 27, which used QCISD/MG3S for CF2 and QCISD(T)/aug-cc-pVQZ for the rest, with unrestricted (restricted) HF references for triplet (singlet) multiplicities, respectively. ph-AFQMC results utilizing CASSCF trial wavefunctions are plotted in Fig. 1, relative to experimental reference values, along with data provided in Ref. 27 for broken-symmetry DFT with the BLYP functional (UKS/BLYP), its spin-projected variant (WA-KS/BLYP), and a composite method which has been shown to produce comparable accuracy to CCSD(T)/CBS69 (W2X). We note that we use a different reference value for NH${}_{2}^{+}$, which is from a genuine experimental measurement70, following Ref. 71.

The ph-AFQMC ST gaps have been converged with respect to active space sizes of the CASSCF trial wavefunctions, an approach which has been shown to produce very high accuracy even for strongly correlated systems.62, 61 However the generation of such trial wavefunctions is, in practice, limited to moderate system sizes due to the procedure’s exponential scaling. In light of applications to large systems such as those found in organic electronics which we will focus on in future investigations, we explore the simplest scalable alternatives, namely UHF and RHF wavefunctions, and single determinants constructed from UB3LYP orbitals. $\langle S^{2}\rangle$ values with respect to the trial wavefunctions are shown in Table 1.

Table 2 provides a rudimentary statistical representation of the accuracy of selected theoretical approaches with respect to experiment, comparing the mean signed error (MSE), mean absolute error (MAE), and maximum error (MaxE). As expected, utilizing CASSCF trial wavefunctions can obtain excellent accuracy, with an MAE of less than a kcal/mol and MaxE of 1.7(3) kcal/mol for O2. For this molecule, we have tried active spaces of 12e8o, 8e12o, and 10e15o. The latter two active spaces produce statistically indistinguishable ST gaps in the CBS limit. Additionally, we verified that the total energies of both the singlet and triplet states in the aug-cc-pvtz basis differ by less than 1$mHa$ going from one active space to the next. This suggests that the AFQMC/CAS result is converged with respect to active space size, and the remaining deviation from experiment is likely due to inaccuracies of the optimized geometries and/or the zero-point energies used to correct the experimental result.

The performance of UKS/BLYP is unsurprisingly poor, given the high level of spin-contamination in the singlet states revealed in Table 1. The notable onset of larger errors in Fig. 1, i.e. for NH, NF, O2, OH+, O, C, and Si, is found to roughly correlate with the presence of spin-contamination in the unrestricted wavefunctions. The Yamaguchi correction clearly improves upon the UKS/BLYP results, however the MAE of 7.3 kcal/mol and MaxE of 17.4 are still very large. W2X is relatively more robust, with an MAE of 3.1 kcal/mol; however, the MaxE of 6.9 kcal/mol illustrates the difficulty in describing biradical systems even with “gold-standard” single-reference methods.

In contrast, the AFQMC/U approach shows a significant improvement in accuracy (all ph-AFQMC results using UHF, UB3LYP, and RHF trials are shown in the Supporting Information). We note that the MaxE for AFQMC/UHF is for the H2CC molecule. The UHF wavefunction for the triplet state, which should largely be of single-reference nature, still exhibits significant spin-contamination, as seen in Table 1. The Slater determinant derived from UB3LYP orbitals appears to be relatively uncontaminated, with an $\langle S^{2}\rangle$ value of 2.03 while still benefiting from the additional variational freedom due to the use of unrestricted orbitals. When this is used as the trial wavefunction for ph-AFQMC, the resulting ST gap in the CBS limit is -49.9(9) kcal/mol, which significantly reduces the deviation from experiment from 3.1(7) to -1.3(9) kcal/mol. While more data points are needed to validate a more general claim, this case suggests that when a single-reference spin-state exhibits spin-contamination, AFQMC/UB3LYP can improve the accuracy of ST predictions over AFQMC/UHF. This is reflected in the protocol specified earlier for AFQMC/U, which utilizes UHF trial wavefunctions for triplet states when $1.9\leq\langle S^{2}\rangle\leq 2.1$, and UB3LYP otherwise.

To make a broader comparison regarding the ability of a variety of electronic structure methods (with similar computational scaling with respect to system size) to predict ST gaps in biradicals, we include calculated values from methods highlighted in Ref. 14 for a subset of 8 biradicals. We plot this data in Fig. 2, and provide a comparative statistical summary in Table 3.

For this subset of cases, there is no distinction between the AFQMC/UHF and AFQMC/U procedures, and AFQMC/U and AFQMC/CAS yield equivalent accuracy, considering statistical error bars. Both produce MAEs of $\sim$1 kcal/mol and MaxEs of $\sim$2 kcal/mol, comparing favorably to all other methods. As shown in the Supporting Information, the accuracy of AFQMC/UB3LYP and AFQMC/RHF is similar to that obtained via spin-flip methods for these systems.

3.2 Benzyne Isomers

In this section we consider the ortho-, meta-, and para-benzyne isomers, shown in Fig. 3, and compare predicted ST gaps with precise, gas-phase experimental measurements. The ground state for all isomers is a singlet, and biradical character correlates with the distance between the unpaired electrons (ortho $<$ meta $<$ para).13. These systems are of scientific interest in their own right, e.g. singlet para-benzyne is a biradical that can abstract hydrogen atoms from specific positions in DNA, potentially enabling antitumor activity.73, 74.

Ab initio calculations are difficult due to the strongly correlated biradical electrons. In particular, singlet para-benzyne exhibits orbital instabilities at the RHF level, and subsequent RHF-based correlation methods produce poor results. Better results are obtained with the broken-symmetry UHF reference, with complete spin symmetry-restoration via, e.g., subsequent CCSD calculation.74 With this in mind, we explore UHF and UB3LYP trial wavefunctions for ph-AFQMC, in addition to multi-determinant trial wavefunctions from CASSCF.

Our ph-AFQMC calculations utilize geometries obtained from SF-DFT/6-311G* using the so-called 50/50 functional13. The ortho- and meta- geometries are provided in Ref. 14, the para- geometry in Ref. 25. In Ref. 14, ST gaps of para-benzyne are based on SF-CCSD/cc-pVTZ optimized geometries25. However we have performed ph-AFQMC calculations on this SF-CCSD geometry with CASSCF trial wavefunctions utilizing increasingly large active spaces (up to 8 electrons in 16 orbitals), and still find a residual error from experiment of 2.9(7) kcal/mol in the cc-pVTZ basis (which is altered by only 1 kcal/mol in the estimated CBS limit). For this reason, and for consistency with the ortho- and meta- isomers, we use the SF-DFT/6-311G* geometry for the para-isomer as well.

It is noteworthy that Refs. 14 and 23 do not account for ZPE contributions in comparing calculated electronic ST energies to experimentally measured quantities. For meta-benzyne this correction is $\sim$1 kcal/mol, so we choose here to subtract out ZPE values, taken from Ref. 13, from the experimental data when comparing with purely electronic predictions.

$\langle S^{2}\rangle$ values of the candidate unrestricted trial wavefunctions are shown in Table 4 for the benzyne isomers. The UHF singlet and triplet states are both significantly contaminated in all isomers, with the singlet state of para-benzyne is severe case ($S^{2}$=1.68). UB3LYP reduces the amount of spin-contamination in all cases, however it does not always eliminate it, e.g. $\langle S^{2}\rangle$ of singlet para-benzyne is reduced to 0.92. The AFQMC/U method will use UB3LYP trial wavefunctions for all isomers, since $\langle S^{2}\rangle_{singlet}^{UHF}>1.1$ for ortho- and para- benzynes, and $\langle S^{2}\rangle_{triplet}^{UHF}>2.1$ for all.

The resulting ST gaps are shown in Table 5, and the deviations of the predicted values from ZPE-corrected experimental results are shown in Fig. 4. AFQMC/UHF and AFQMC/UB3LYP values are shown separately in the Supporting Information.

SF-oo-CCD and SF-CCSD(T) both consistently obtain sub kcal/mol accuracy, however they are computationally infeasible for larger systems due to the respective $O(M^{6})$ and $O(M^{7})$ scaling. We thus focus presently on methods with lower scaling.

AFQMC/CAS and AFQMC/U produce predictions of comparable accuracy, with maximum errors just outside 1 kcal/mol. AFQMC/UHF achieves good accuracy for the ortho- and meta- isomers, however the severely spin-contaminated singlet state of para-benzyne results in a relatively large overestimation of the ST gap. For this latter system AFQMC/UB3LYP substantially reduces the deviation from experiment from 5.7(8) to 1(1) kcal/mol, and we note that similarly pronounced corrections are observed in pp-RPA calculations when the B3LYP reference for the ($N-2$)-electron system is used instead of HF.14 AFQMC/U and AFQMC/UB3LYP are equivelent for this set of molecules, and thus we again find that AFQMC/U, utilizing only unrestricted single-determinant trials, can produce results of comparable accuracy to both AFQMC/CAS and experiment.

We must emphasize the importance of explicitly breaking the spin-symmetry when obtaining UHF and UB3LYP trial wavefunctions, and choosing the solution with lowest-energy. For instance, the calculated ST gap using a spin-pure ($\langle S^{2}\rangle$ = 0) B3LYP trial for para-benzyne is -15 kcal/mol, vs +2.2 kcal/mol when the lower-energy $\langle S^{2}\rangle=0.92$ unrestricted singlet state is used.

As a final remark in this section, we note that while indeed the errors from UB3LYP (i.e. without subsequent QMC) are generally reduced in comparison with those from the small molecule biradical set, we still find significant errors of -8.7, -5.8, and -0.9 kcal/mol for the ortho-, meta-, and para-benzynes. UB3LYP systematically underestimates the ST gaps, due to the unrealistically low energy of the broken-symmetry singlet state. Rather unexpectedly, however, the magnitude of the errors here are inversely correlated with diradical character.

These results suggest that for the benzyne isomers, which exhibit strong biradical character while sharing features similar to the planar aromatic ring systems relevent to chemical photocatalysis, the ST gaps are accurately predicted by ph-AFQMC with both CASSCF and single-determinant trials, and in the cc-pVTZ basis. The accuracy of AFQMC/UHF, even in its spin projected form, is compromised by the heavily spin-contaminated singlet state in para-benzyne, though AFQMC/UB3LYP provides an improved prediction, and is utilized in our AFQMC/U method.

3.3 Polyacenes

Having shown that AFQMC/U can accurately describe molecules with strong biradical nature, we now show that this computational approach can scale to larger molecules, focusing on polyacenes from naphthalene to pentacene, shown in Fig. 3. These molecules all have singlet ground states. We use geometries from Ref. 77 for the acenes, which were computed at the unrestricted B3LYP/6-31G(d) level of theory.

Table 6 shows that the extent of spin-contamination in the singlet UHF states increases with the number of fused rings. In contrast, wavefunctions constructed from UB3LYP orbitals are spin-pure, consistent with previous computational studies78, 33, 32.

ST gaps calculated with ph-AFQMC in the cc-pVTZ basis are shown in Table 7, alongside predictions from state-of-the-art ab initio methods and available experimental data. When more than one experimental value is given in Ref. 66, we choose the value that is closest to that shown in Ref. 42. Zero point energy corrections, which are subtracted out of the experimental values, are required in order to compare calculated electronic energies with experiment, and we utilize numbers from Ref. 66 derived from B3LYP/6-31G(d) geometry optimizations and frequency calculations.

It is important to recognize that all the experimental values, except in the case of anthracene, cannot be fairly compared directly with gas-phase calculations. Ref. 33 conveniently provides details of many of the experimental measurements, which are reproduced here. While the adiabatic ST gap of anthracene was obtained via gas-phase photoelectron spectroscopy, the reported experimental measurement for naphthalene was done in ether-isopentane-alcohol (solid) solvents at 77K. The tetracene measurement was done in poly(methyl methacrylate) matrix at 298K, and pentacene was measured in a tetracene matrix at 298K. Clearly, the gas-phase 0K conditions assumed in our calculations are not consistent with the realistic experimental conditions for most of the polyacenes studied here.

We observe that the ph-AFQMC predictions are insensitive to the trial wavefunction used for these polyacene systems, and show results from UHF, RHF, and UB3LYP trial wavefunctions in the Supporting Information. Given that naphthalene through pentacene exhibit little biradical character42, the extreme spin-contamination of the UHF solutions shown in Table 6 is likely not representative of strong electron correlation, and hence restricted trial wavefunctions can also be expected to yield accurate results. Moreover, all ph-AFQMC variants unambiguously achieve high accuracy with respect to the gas-phase measurement for anthracene. The ph-AFQMC predictions for the other polyacenes, which were experimentally probed in (solid) solvent matrix, appear to systematically overestimate the experimental values by a few kcal/mol. For naphthalene, we performed ph-AFQMC calculations with CASSCF(10e,10o) trial wavefunctions, which gave a ST gap of 67(1) kcal/mol. This value is in agreement with that from AFQMC/U, given statistical error bars, and lies above the reported ZPE-corrected experiment by some 3 kcal/mol. Very similar overestimations by AFQMC/U are seen for $n$=2,4,5 and are corroborated by CCSD(T)/FPA, B3LYP/pp-RPA, and DMRG-pDFT (though not by UB3LYP and ACI-DSRG-MRPT2). Thus it seems reasonable to hypothesize that the calculations’ neglect of molecular environment may be responsible, though admittedly there are a number of other factors that might be expected to contribute. For instance, adiabatic ST gaps are known to be sensitive to the optimized geometries, which can result in variations on the order of 1-3 kcal/mol.66

Overall, given the uncertainties due to the treatment of temperature, solvent, and molecular geometries in these acene calculations, ph-AFQMC with single-determinant trial wavefunctions gives satisfactory agreement with both experiments and other highly-accurate electronic structure methods, and moreover can scale with near-perfect parallel efficiency to systems as large as pentacene in a triple-zeta basis, which has 146 electrons and 856 basis functions.

4 Conclusions

With CASSCF trial wavefunctions, ph-AFQMC can predict ST gaps with sub kcal/mol accuracy with respect to gas-phase experimental measurements for a set of 13 small molecules with singlet states of strong biradical character. However, for large systems the generation of such a trial wavefunction quickly becomes impractical. The main result of this work is that near-chemical accuracy for gas-phase ST gaps can also be obtained, even for strongly correlated biradical systems, with a spin-projected ph-AFQMC technique, which initializes walkers with restricted determinants while using an unrestricted single-determinant trial wavefunction to implement the phaseless constraint. We establish a quantitative criteria for choosing UHF or UB3LYP orbitals based on the spin-contamination of the UHF wavefunction, and the resulting AFQMC/U methodology is validated on the small molecule test set, the ortho-, meta-, and para-benzyne isomers, and of all the polyacenes for which experimental results are reported (though a true gas-phase experiment is only available for anthracene). The ph-AFQMC method is shown to provide a balanced and robust approach to accurately predicting ST gaps for molecules relevant to chemical photocatalysis.

With computational cost scaling as $O(M^{4})$ or lower,79 and with near-perfect parallel efficiency,62, 61 these calculations take only minutes of wall-time on GPU-accelerated supercomputing resources. While further investigation will be required to realistically treat solvation effects, this study paves the way for future work that will focus on a large set of both known and potential photocatalysts. Another possible extension we plan to undertake in the future is the computation of spin gaps in strongly correlated solids.

{acknowledgement}

J.S. acknowledges Roel Tempelaar, Joonho Lee, Mario Motta, Hao Shi, Ian Dunn, Andrew Pun, and Zack Strater for insightful discussions. This research used resources of the Center for Functional Nanomaterials, which is a U.S. DOE Office of Science Facility, and the Scientific Data and Computing Center, a component of the Computational Science Initiative, at Brookhaven National Laboratory under Contract No. DE-SC0012704. An award of computer time was provided by the INCITE program. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. The Flatiron Institute is a division of the Simons Foundation. D.R.R. acknowledges support from NSF CHE 1839464, and S.Z. from DOE DE-SC0001303.

{suppinfo}