Fixed node-diffusion Monte Carlo achieves chemical accuracy in predicting substituent effects on activation energies and reaction enthalpies for methyl radical addition to substituted olefins

Abstract

Computing accurate activation barriers and reaction enthalpies is essential for the development of kinetic mechanisms and prediction of reaction outcomes. However, the computationally intensive nature of accurate quantum calculations and lack of experimental data present challenges. In this study, eighteen radical addition reactions relevant to free radical polymerization were used to assess the accuracy of single-determinant fixed-node diffusion Monte Carlo (FN-DMC) in predicting activation barriers and reaction enthalpies. Using CCSD(T)/aug-cc-pVTZ as a reference, FN-DMC acquired a mean absolute deviation (MAD) of 4.5(5) kJ mol⁻¹ for activation energies and 3.3(5) kJ mol⁻¹ for reaction enthalpies. Activation energies and reaction enthalpies relative to the unsubstituted olefin exhibited relative MAD values of 1 kJ mol⁻¹ and accurately predicted trends in substituent effects within statistical margins, except for the activation energies of the butanenitrile and 2-methoxybutyl radical systems and the reaction enthalpies of the ethoxypropyl and butyl analogs. With the demonstrated accuracy, favorable computational scaling, and highly parallelizable nature of FN-DMC, it should be feasible to use FN-DMC to investigate activation energies and reaction enthalpies for larger systems, such as oligomers, where coupled-cluster based methods may be limited by computational resources.

---

1 Introduction

Free radical polymerization (FRP) plays an important role in the global production of synthetic polymers, accounting for approximately half of the total output.^1,2 Due to its versatility and high reactivity, conventional FRP techniques have been used to manufacture a variety of commodities such as polyethylene, polystyrene, polyvinylacetate, and polybutadienes.³ In its simplest form, FRP is made up of initiation, propagation, chain transfer and termination reactions. The polymer elongates its chain by forming new sigma bonds through an attack of a free radical on a pi bond. Consequently, the resulting polymer retains a radical within its structure, allowing the process to repeat with excess monomer. A primary limitation of FRP is the extensive chain termination and chain transfer processes, leading to broad molecular weight distributions.² Advances in controlled polymerization techniques, such as atom transfer radical polymerization (ATRP) and reversible-deactivation radical polymerization (RDRP), have significantly improved structural control by mitigating radical concentrations.⁴ However, to unlock their full potential for intricate polymer architectures, precise kinetic data is essential. These processes may involve up to tens of kinetically distinct chemical reactions, each contributing to the overall rate of polymerization and molecular weight.⁵ These propagation rate constants are typically investigated using pulsed-laser polymerization in conjunction with size exclusion chromatography (PLP-SEC), as recommended by IUPAC.⁶ By introducing short laser pulses onto a monomer solution to initiate and control the polymerization process, precise control over the reaction conditions and resulting polymer properties are enabled. Rate constants have been determined for a wide range of monomers in various solvents to elucidate structure–reactivity relationships.^7,8 Despite its utility, the technique faces limitations with slowly terminating monomers and chain transfer reactions resulting in weak signals and broad molecular weight distributions.⁹ In turn, these diffusion-controlled reactions are further studied using single-pulsed laser polymerization electron paramagnetic resonance spectroscopy (SP-PLP-EPR).⁹

With rapid advancements in computing power and ab initio quantum chemical methods, first-principles predictions have become a vital precursor to predicting kinetic parameters within chemical accuracy (1 kcal mol⁻¹ = 4.184 kJ mol⁻¹). With the assistance of transition state theory,¹⁰ these quantum chemistry methods are able to predict molecular geometries, transition states, and activation energies. The “gold standard” for gas-phase systems well described by a single-determinant wavefunction is coupled cluster with single, double, and perturbative triple excitations [CCSD(T)].^11,12 It is typically only applied to the smallest systems due to its substantial computational scaling of O(N⁷), where N is a measure of the system's size (the number of atoms, assuming a fixed number of basis functions). One strategy to overcome this limitation is to use a proxy for canonical CCSD(T), such as the domain-based local pair natural orbital coupled cluster [DLPNO-CCSD(T)] approach.^13–16 This technique employs localized occupied orbitals and includes only electron pairs that significantly contribute to correlation energy, while the remaining pairs are treated at the MP2 level. Its performance relative to its canonical counterpart for radical addition reactions have been evaluated within numerous studies.^17–20 Recently, we assessed the performance of DLPNO/CCSD(T) using both Hartree–Fock and Kohn–Sham reference orbitals for a set of 18 radical addition reactions.²¹ The findings indicated that DLPNO-CCSD(T)/CBS with unrestricted Hartree–Fock reference orbitals and NormalPNO parameters accurately computed activation energies and reaction enthalpies, achieving mean absolute deviations of 3.5 kJ mol⁻¹ and 1.4 kJ mol⁻¹, respectively, when compared to CCSD(T)/aug-cc-pVTZ.

Alternatively, quantum composite methods are able to obtain results comparable to higher levels of theory by combining larger basis sets with modest levels of theory and smaller basis sets with more demanding theory. The complete basis set methods (CBS),^22–24 Weizmann-n theory,^25–27 Guassian-n theory,^28–31 and correlation consistent Composite Approaches (ccCA)^32–34 have been routinely used and have demonstrated chemical accuracy for thermodynamic quantities.^35–40 Although these methods have a reduced computational overhead compared to CCSD(T) due to the use of smaller basis sets, they remain impractical for large systems, such as chain transfer to oligomers or trimers involving more than forty atoms.

Density functional theory (DFT) has long been favored as an alternative to wavefunction theory due to its ability to account for electron correlation across a broad spectrum of properties and systems at a relatively low computational cost. However, the accuracy of DFT depends on the chosen exchange–correlation approximation.⁴¹ A density functional approximation (DFA) that performs well for one system may not be effective for other systems. Recently, Duan and colleagues developed a “transferable recommender” system to select the most suitable DFA for specific systems and properties.⁴² The model, when applied to vertical spin splitting energies, achieved an accuracy of 2 kcal mol⁻¹. Such recommender systems have yet to be extensively applied to kinetics and thermochemistry. Despite advancements in DFT, accurately computing thermochemical and kinetic properties remains challenging for many functionals due to non-systematic errors.⁴³ Bond dissociation energies, in particular, have proven difficult for various DFAs.^44–46 A comprehensive database, MDCDDB84, which contains reaction enthalpies and barrier heights, demonstrated that most hybrid and double-hybrid density functionals fail to achieve chemical accuracy, with the exception of ωb97M(2).⁴⁷ Radical addition reactions involving delocalization of the unpaired electron, such as those with nitrile or benzyl functional groups, show discrepancies greater than 10 kJ mol⁻¹ compared to experimental values.^48,49 These open-shell species present additional challenges due to their multi-determinant character and spin contamination. Multiconfiguration pair-density functional theory (MC-PDFT) addresses these issues by combining wavefunction theory for static correlation and DFT for dynamic correlation. In MC-PDFT, the Coulombic and electronic kinetic energies are computed from a reference multiconfiguration wavefunction, while the remaining energy components are derived from a density functional, referred to as the on-top density functional.^50,51 Numerous studies have explored the application of MC-PDFT to activation energies, highlighting its potential to improve the accuracy of DFT for complex systems.^52–56

With the increase of computational power and access to large clusters, fixed-node diffusion Monte Carlo has emerged as a promising alternative to these deterministic methods. FN-DMC, a stochastic member of the quantum Monte Carlo (QMC) family, represents the wavefunction as an ensemble of electron-position configurations, termed walkers, sampled from an importance-weighted probability distribution.^57,58 The fundamental principle involves imposing an imaginary time exponential projection to extract the ground state from an initial trial wavefunction. Proper discretization of this propagation is crucial to accurately simulate the continuous evolution described by the imaginary time Schödinger equation. Recently, a novel technique has been developed to enable the utilization of larger timesteps, alleviating the necessity for extrapolating to infinitesimally small timesteps.⁵⁹ To accomodate the fermionic nature of electrons, the nodal surface of the trial wavefunction remains fixed, thereby preventing walkers from crossing these nodes. Using simple Slater–Jastrow trial wavefunctions, DMC exhibited exceptional accuracy for small weakly bound systems.⁶⁰ Its scalability has led to successful simulations of noncovalent bonds in large supramolecular complexes containing hundreds of atoms, such as the DNA-intercalating anticancer drug ellipticine complex⁶¹ and the buckyball catcher⁶² from the S12L database.⁶³ However, when compared to CCSD(T) for polarizable supramolecules, discrepancies of 8 kcal mol⁻¹ were observed.⁶⁴ Employing single-determinant nodal surfaces for covalent species, FN-DMC exhibited good agreement with experimental barrier heights for both hydrogen transfer reactions⁶⁵ and non-hydrogen transfer reactions,⁶⁶ with mean absolute deviations of 1.0(3) kcal mol⁻¹ and 1.5(1) kcal mol⁻¹, respectively. Krongchon et al. identified a correlation between small HOMO–LUMO gaps and the magnitude of fixed-node error,⁶⁶ consistent with previous findings emphasizing the necessity of multi-determinant character for transition states with stretched bonds.⁶⁷ Similar trends were observed for open-shell bond dissociation energies, which exhibited large discrepancies,⁶⁸ while closed-shell precursor bond dissociation energies were modeled effectively.^68,69 While multi-determinant trial wavefunctions have yielded improvements in atomization energies within the G2 test set,^70,71 single-determinant wavefunctions have nonetheless provided chemical accuracy for a broad range of properties, including radical stabilization energies,⁷² neutral electron affinities and ionization potentials of second and third row atoms and molecules (with ionic species showing significant improvement with a multi-determinant wavefunction),⁷³ and hydrogen abstraction reaction enthalpies.²¹ The current study is an extension of our previous work, which compared FN-DMC-derived reaction enthalpies for radical addition and hydrogen abstraction reactions with those obtained using composite quantum chemical methods.⁷⁴ Our present work aims to further evaluate the capabilities and limitations of FN-DMC by applying it to calculate reaction enthalpies and activation energies of eighteen radical addition reactions to substituted olefins and comparing the results to coupled cluster calculations for the same reactions.²¹

2 Computational details

To assess the performance of FN-DMC in generating activation energies and reaction enthalpies pertinent to FRP, eighteen radical addition reactions to substituted olefins from Radom and Fischer's comprehensive test set⁷⁵ are examined. Sixteen of these reactions involve methyl radical addition, while the remaining two involve hydroxymethyl (MOH) and cyanomethyl (MCN) radical additions to ethylene. Experimental activation energy data for the phenyl-substituted olefin is unavailable, thus only the reaction enthalpy is considered. Substituent abbreviations are as follows: Me = methyl, Et = ethyl, OEt = ethoxy, OAc = acetoxy, and Ph = phenyl.

The optimized geometries are the same as those reported in ref. 74 obtained using B3LYP^76–79/cc-pVTZ.⁸⁰ For comparative analysis, we include experimental values, G3(MP2)-RAD, and CBS-RAD from ref. 75, and CCSD(T)/aug-cc-pVTZ and DLPNO-CCSD(T)/CBS from ref. 21.

The computational protocol followed here is similar to that previously tested and optimized for methyl radical addition to ethylene.⁷⁴ Specifically, the FN-DMC trial wavefunctions used for importance sampling is of the Slater–Jastrow form,

ψ_T(x⃑) = ψ_A(x⃑)e^J(x⃑) (1)

where ψ_A(x⃑) is an antisymmetric Slater determinant from Hartree–Fock or density functional theory, and J(x⃑) is a Jastrow factor containing explicit electron correlation terms. This work includes terms for electron–nucleus, electron–electron, and electron–electron–nucleus interactions. The Jastrow factors are optimized using the linear method⁸¹ with a cost function comprising a 95 : 5 mixture of re-weighted variance and energy. Orbital optimization is not performed, as its unfavorable scaling with basis set size did not justify the minor improvements observed for olefin-related systems using pseudopotentials.^71,82 Using PySCF,^83–85 a single-determinant trial wavefunction is constructed using PBE0/BFD-VTZ single-particle orbitals alongside the Hartree–Fock energy-consistent pseudopotential, Burkatzki–Filippi–Dolg (BFD).⁸⁶ Throughout this work, FN-DMC refers to simulations using this PBE0/BFD-VTZ trial wavefunction. While multi-determinant wavefunctions can improve FN-DMC accuracy,^71,87 previously published work indicates that the systems examined here are predominantly single-determinant in nature.²¹ This conclusion is supported by coupled cluster T₁ diagnostics,⁸⁸ despite some spin contamination at the UHF level in transition states. Casula's v1 size-consistent T-moves approach is employed to handle the nonlocal pseudopotential beyond the locality approximation.⁸⁹ DMC calculations are carried out using 2048 walkers and a timestep of 0.025 a.u. with a Green's function propagator. This timestep is chosen based on prior convergence tests with ECPs, which showed deviations within 1 mHa of the extrapolated zero-timestep limit.⁷⁴ Energy differences further reduce residual timestep bias to statistically insignificant levels. The potential walker population bias is evaluated using smaller walker counts, confirming that 2048 walkers were sufficient to avoid such bias.⁷⁴ After 4000 equilibrium steps, DMC statistics are accumulated using a sufficient number of Monte Carlo steps, ensuring that statistical errors in activation energies and reaction enthalpies were less than 1.0 kJ mol⁻¹. All VMC and DMC calculations are carried out using QMCPACK code v3.12.0⁸² on Northern Illinois University's Gaea hybrid compute cluster.

3 Results

3.1 Activation energies

The activation energies computed using FN-DMC are summarized in Table 1. Compared to the experimental values obtained from kinetic ESR measurements,^75,90 FN-DMC exhibited a mean absolute deviation (MAD) of 5.8(5) kJ mol⁻¹. This value is in close agreement with those obtained from DLPNO-CCSD(T)/CBS and G3(MP2)-RAD, differing by less than 1 kJ mol⁻¹. In contrast, the canonical CCSD(T)/aug-cc-pVTZ and CBS-RAD methods yielded MAD values of 2.9 kJ mol⁻¹ and 2.2 kJ mol⁻¹ respectively, suggesting improved accuracy. Given the substantial experimental error of ±5 kJ mol⁻¹, this comparison should be interpreted with caution. Consequently, the subsequent analyses adopt CCSD(T)/aug-cc-pVTZ with UB3LYP reference orbitals²¹ as the reference method. With this new reference, the MAD of FN-DMC was 4.5(5) kJ mol⁻¹, which is comparable to the performance of the other post-Hartree–Fock methods. Although CBS-RAD exhibited the smallest MAD, it is important to note that this statistic was calculated with 5 fewer data points, including the particular challenging 2-methylbutanenitrile radical (X = Me, Y = CN). In addition, FN-DMC consistently overestimates the activation energies across the dataset, highlighting a systematic trend. To aid in the interpretation of these results, a box-and-whisker plot of absolute residual errors relative to CCSD(T)/aug-cc-pVTZ is presented in Fig. 1a. The median absolute residuals (horizontal line) for all methods are within ±2 kJ mol⁻¹, with FN-DMC showing the most consistent results, as evidenced by the narrowest interquartile range (box height) and shortest whiskers. FN-DMC exhibited one significant outlier, 2-methylbutanenitrile represented by an individual data point. In comparison, DLPNO-CCSD(T)/CBS displays the lowest central tendency but exhibits a positive skew in its residuals.

Table 1 A comparison of experimental and calculated activation energies (kJ mol⁻¹) at 0 K for the addition of a methyl radical to olefins (CH₂ = CXY). The experimental values are reported with a ±5 kJ mol⁻¹ error bar. The statistical uncertainty of the last digit for FN-DMC is indicated in parentheses. Statistical analysis was performed using both CCSD(T)/aug-cc-pVTZ and experimental values (in brackets) as references

X^a	Y^a	Exp.^b	CCSD(T)/aug-cc-pVTZc	FN-DMC		DLPNO/CBS^c		G3(MP2)-RAD^b		CBS-RAD^b
a Substituent abbreviations: Me = methyl, Et = ethyl, OEt = ethoxy, OAc = acetoxy. b Values taken from ref. 75 and 90. c Values taken from ref. 21. d Activation energies for the CH2OH addition to ethylene. e Activation energies for the CH2CN addition to ethylene. f Mean deviation (MD), mean absolute deviation (MAD), maximum absolute deviation (MaxAD), standard deviation (SD). g Relative values were calculated as the difference between activation energy for the given reaction and the corresponding value for the unsubstituted olefin (X = H, Y = H).
H	H	30.4	32.9	36.6(3)		35.7		36.2		29.7
H	Me	29.0	31.8	35.0(5)		33.9		34.8		27.6
H	Et	27.5	31.4	34.9(5)		33.6		34.5		26.8
Me	Me	26.8	29.7	33.8(5)		31.8		32.8		25.3
Me	OMe	28.2	30.9	37.0(5)		36.8		36.2		—
H	OEt	28.3	37.2	41.8(5)		39.9		34.8		—
H	OAc	29.4	29.3	33.6(5)		32.7		32.4		—
H	Cl	25.2	26.1	29.8(4)		28.4		29.2		22.2
H	SiH3	24.8	26.8	30.4(4)		29.1		29.9		23.1
Me	Cl	23.1	25.1	28.7(5)		27.5		27.8		20.4
Cl	Cl	17.4	19.3	25.2(4)		22.0		22.2		16.7
H	CO2Me	17.6	14.5	21.0(5)		21.1		21.1		—
H	CN	16.4	16.1	20.0(4)		20.2		21.1		13.8
H	CHO	16.0	11.7	16.1(5)		16.0		22.6		14.6
Me	CN	15.3	15.5	22.4(5)		20.1		20.7		—
MOH^d		36.9	28.3	33.2(4)		32.5		34.1		27.6
MCN^e		30.0	31.4	34.9(4)		36.4		35.7		28.8
Absolute values^f
MD		—	[0.9]	4.5(5)	[5.4(5)]	3.5	[4.4]	4.2	[5.1]	−2.8	[−2.2]
MAD		—	[2.9]	4.5(5)	[5.8(5)]	3.5	[4.9]	4.2	[5.5]	3.3	[2.2]
MaxAD		—	[8.9]	6.9(5)	[13.5(5)]	6.6	[11.6]	10.9	[10.0]	4.7	[9.3]
SD		—	[3.8]	1.2	[3.7]	1.4	[3.4]	2.2	[2.7]	2.1	[2.4]
Relative values^f^g
MD		—	[−1.6]	0.9(5)	[−0.7(5)]	0.9	[−0.8]	1.1	[−0.6]	0.4	[−1.5]
MAD		—	[2.7]	1.0(5)	[2.5(5)]	1.2	[2.2]	1.6	[1.8]	1.4	[1.5]
MaxAD		—	[11.2]	3.3(5)	[9.9(4)]	3.8	[9.7]	7.6	[8.6]	6.1	[8.6]
SD		—	[3.8]	1.2	[3.7]	1.4	[3.4]	2.2	[2.7]	2.1	[2.4]

---

Fig. 1 (a) A box-and-whisker plot illustrating residual errors, calculated as the difference between predicted and CCSD(T)/aug-cc-pVTZ reference activation energies. The box height represents the interquartile range (IQR), encompassing the middle 50% of errors, with the horizontal line indicating the median. Whiskers extend to the full data range excluding outliers, which are shown as circles beyond 1.5 times the IQR. This visualization highlights how well each method predicts activation energies and identifies systematic trends and inconsistencies. (b) A relative ranking matrix of activation energies, calculated as the difference between each reaction and the unsubstituted olefin, showcasing the ability of post-Hartree–Fock methods to capture substituent effects. CCSD(T)/aug-cc-pVTZ values (the first column) determine reference colors, with darker colors and higher ranks corresponding to larger relative activation energies (e.g. substituents leading to a larger activation barrier). Colors in subsequent columns differing from those in the first column give a vivid visual indication of rankings differing from those determined by CCSD(T)/aug-cc-pVTZ. White spaces indicate that no results are available.

For chemists, the ability for a computational method to predict relative substituent effects on activation energies is of primary importance. This capability can be assessed by comparing the activation energy of a given reaction with that of the unsubstituted olefin (X = H, Y = H). As shown in Table 1, FN-DMC achieved a MAD of 1.0(5) kJ mol⁻¹ for the relative activation energies, with the 2-methylbutanenitrile radical (X = Me, Y = CN) exhibiting the highest absolute maximum deviation (MaxAD) of 3.3(5) kJ mol⁻¹. A relative ranking matrix (Fig. 1b) highlights the ranking of relative differences in activation energies, with higher ranks (depicted by darker cells and larger numerical values) representing greater relative differences, making outliers easily identifiable. Among the tested methods, FN-DMC demonstrated the closest agreement with CCSD(T)/aug-cc-pVTZ, misranking only the butanenitrile (X = H, Y = CN) and 2-methoxybutyl radical (X = Me, Y = OMe) systems.

3.2 Reaction enthalpies

In addition to understanding the barrier heights associated with radical polymerization reactions, evaluating the thermodynamic feasibility is equally important. Consequently, the reaction enthalpies were computed with FN-DMC and compared with the results of four post-Hartree–Fock methods (Table 2). Similar to the activation energy comparisons, caution must be advised when comparing reaction enthalpies with experimental values due to a quoted experimental error of ±10 kJ mol⁻¹. Using CCSD(T)/aug-cc-pVTZ as the reference, FN-DMC was shown to predict chemically accurate reaction enthalpies, achieving a MAD of 3.3(5) kJ mol⁻¹. Fig. 2a presents the residual errors of FN-DMC and three post-Hartree–Fock methods relative to CCSD(T)/aug-cc-pVTZ in a box-and-whisker plot. Notably, DLPNO-CCSD(T)/CBS and CBS-RAD exhibit better performance, with their median (Q2 quartile) values closer to zero. FN-DMC demonstrates remarkable consistency, as evidenced by its small interquartile range, while G3(MP2)-RAD underperforms, showing large residual errors. The 1,1-dichloropropyl radical system (X = Cl, Y = Cl) proves particular challenging for FN-DMC, exhibiting a MaxAD of 5.2(5) kJ mol⁻¹. Notably, FN-DMC overestimates the activation energy for the 2-methoxybutyl radical by 6.1(5) kJ mol⁻¹, while accurately predicting the reaction enthalpy, with an overestimation of only 0.4(5) kJ mol⁻¹.

Table 2 A comparison of experimental and calculated reaction enthalpies (kJ mol⁻¹) at 0 K for the addition of a methyl radical to olefins (CH₂ = CXY). The experimental values are reported with a ±10 kJ mol⁻¹ error bar. The statistical uncertainty of the last digit for FN-DMC is indicated in parentheses. Statistical analysis was performed using both CCSD(T)/aug-cc-pVTZ and experimental values (in brackets) as references

X^a	Y^a	Exp.^b	CCSD(T)/aug-cc-pVTZc	FN-DMC		DLPNO/CBS^c		G3(MP2)-RAD^b		CBS-RAD^b
a Substituent abbreviations: Me = methyl, Et = ethyl, OEt = ethoxy, OAc = acetoxy, Ph = phenyl. b Values taken from ref. 75. c Values taken from ref. 21. d Reaction enthalpies for the CH2OH addition to ethylene. e Reaction enthalpies for the CH2CN addition to ethylene. f Computed with cc-pVTZ basis set.^21 g Mean deviation (MD), mean absolute deviaion (MAD), maximum absolute deviation (MaxAD), standard deviation (SD). h Relative values were calculated as the difference between reaction enthalpy for the given reaction and the corresponding value for the unsubstituted olefin (X = H, Y = H).
H	H	−93.3	−93.2	−96.4(3)		−92.0		−86.3		−91.3
H	Me	−92.2	−94.5	−98.7(4)		−94.8		−86.7		−93.6
H	Et	−94.2	−95.0	−98.6(5)		−95.3		−87.1		−94.6
Me	Me	−93.6	−90.7	−93.6(4)		−91.0		−82.9		−91.6
Me	OMe	−89.2	−87.8	−88.2(5)		−84.0		—		—
H	OEt	−95.1	−95.2	−97.9(5)		−94.7		−86.7		—
H	OAc	−95.9	−113.7	−117.9(6)		−112.4		−92.7		—
H	Cl	−100.8	−107.4	−112.0(4)		−107.3		−98.8		−107.1
H	SiH3	−95.8	−102.8	−106.8(4)		−103.0		−95.0		−100.6
Me	Cl	−93.6	−105.5	−109.7(4)		−105.1		−96.9		−108.3
Cl	Cl	−119.6	−124.2	−129.3(5)		−123.0		−114.1		−124.6
H	CO2Me	−111.5	−119.0	−120.8(5)		−114.6		−107.4		—
H	CN	−133.3	−126.6	−129.4(4)		−125.5		−118.7		−130.4
H	CHO	−116.0	−133.1	−137.5(4)		−131.4		−112.4		−122.9
Me	CN	−125.0	−130.9	−132.1(5)		−127.9		−120.7		—
MOH^d		−81.5	−84.0	−87.3(4)		−81.7		−79.6		−82.6
MCN^e		−65.7	−72.5	−76.5(5)		−69.0		−64.7		−69.1
H	Ph	−144.0	−131.9^f	−134.8(6)		−134.0		—		—
Absolute values^g
MD		—	[−4.7]	−3.3(5)	[−8.0(5)]	1.4	[−2.6]	9.8	[4.8]	1.1	[−3.1]
MAD		—	[5.7]	3.3(5)	[8.1(5)]	1.5	[5.6]	9.8	[5.2]	2.4	[4.2]
MaxAD		—	[17.8]	5.2(5)	[22.0(6)]	4.4	[16.5]	21.0	[14.6]	10.2	[14.7]
SD		—	[6.6]	5.2	[7.2]	4.4	[7.1]	4.6	[4.3]	3.5	[4.9]
Relative values^g^h
MD		—	[−4.9]	−0.2(5)	[−5.1(5)]	0.4	[−4.7]	3.1	[−2.4]	−0.8	[−5.1]
MAD		—	[6.2]	1.2(5)	[6.8(5)]	1.2	[6.1]	3.4	[4.0]	2.5	[5.2]
MaxAD		—	[17.9]	2.9(6)	[17.9(6)]	3.2	[17.8]	14.1	[10.3]	8.3	[16.7]
SD		—	[6.7]	1.5	[7.3]	1.5	[6.7]	4.7	[4.5]	3.6	[5.2]

---

Fig. 2 (a) A box-and-whisker plot illustrating residual errors, calculated as the difference between predicted and CCSD(T)/aug-cc-pVTZ reference reaction enthalpies. The box height represents the interquartile range (IQR), encompassing the middle 50% of errors, with the horizontal line indicating the median. Whiskers extend to the full data range excluding outliers, which are shown as circles beyond 1.5 times the IQR. This visualization highlights how well each method predicts reaction enthalpies and identifies systematic trends and inconsistencies. (b) A relative ranking matrix of relative reaction enthalpies, calculated as the difference between each reaction and the unsubstituted olefin, demonstrating the ability of post-Hartree–Fock methods to capture substituent effects. CCSD(T)/aug-cc-pVTZ values (the first column) determine reference colors, with darker colors and higher ranks corresponding to larger relative reaction enthalpies (e.g. less exothermic substituents). Colors in subsequent columns differing from those in the first column give a vivid visual indication of rankings differing from those determined by CCSD(T)/aug-cc-pVTZ. White spaces indicate that no results are available.

Focusing on the relative substituent effects for reaction enthalpies, Table 2 highlights that FN-DMC achieves exceptional accuracy with a MAD of 1.2(5) kJ mol⁻¹, matching or exceeding the performance of the other post-Hartree–Fock methods. To further evaluate FN-DMC's ability to predict relative substituent effects individually, a relative ranking matrix depicting the ranked relative energy differences of substituted olefins is shown in Fig. 2b. Within this matrix, higher ranks and darker cell colors correspond to greater relative energy differences. FN-DMC correctly ranked the majority of systems (relative to CCSD(T)/aug-cc-pVTZ), with exceptions observed for the ethoxypropyl (X = H, Y = OEt) and butyl (X = H, Y = Me) radical systems. These misrankings can likely be attributed to FN-DMC's statistical error of ±0.5 kJ mol⁻¹, as the energy differences for these systems fall within this range (see SI).

4 Conclusions

In this work, the performance of FN-DMC for predicting activation energies and reaction enthalpies was assessed using a test set of eighteen substituted olefin systems originally proposed by Radom and Fischer.⁷⁵ Given the significant uncertainties in available experimental data, CCSD(T)/aug-cc-pVTZ with B3LYP reference orbitals served as the benchmark. FN-DMC achieved MAD values of 4.5(5) kJ mol⁻¹ for activation energies and 3.3(5) kJ mol⁻¹ for reaction enthalpies. Both thermokinetic properties were consistently overestimated, suggesting a systematic bias. This behavior is likely related to differences in nodal surface quality, as FN-DMC does not benefit from the same degree of error cancellation between species as CCSD(T). Subtle correlation and electron delocalization effects in radicals and transition states may contribute to imbalanced fixed-node errors. Despite these trends, no clear pattern emerged regarding which specific systems exhibited lower FN-DMC accuracy. While a multi-determinant trial wavefunction may improve the nodal surface for the activation energy of the delocalized 2-methylbutanenitrile radical, prior spin-contamination studies²¹ suggest that outliers such as the activation energy of the methyl butyrate radical and reaction enthalpy of the chloropropyl radical are unlikely to be significantly affected due to their strong single-determinant character. Future work examining correlations between nodal surface quality and HOMO–LUMO gaps, as shown in previous work with barrier heights,⁶⁶ may provide further insight. Additional improvements could arise from optimizing orbitals in the presence of the Jastrow factor, which may reduce errors for more delocalized transition states.

Activation energies and reaction enthalpies relative to the unsubstituted olefin, calculated using FN-DMC, exhibited MAD values near 1 kJ mol⁻¹, with uncertainties overlapping this threshold. This stochastic method successfully captured trends in substituent effects within statistical error margins, except for the activation energies of the butanenitrile and 2-methoxybutyl radical systems, and the reaction enthalpies of the ethoxypropyl and butyl analogs.

While DLPNO-CCSD(T) is efficient for small to medium systems, including those studied here, its basis-set incompleteness for larger systems²⁰ will incur significant computational overhead. In contrast, FN-DMC scales as O(N³–N⁴), parallelizes efficiently (especially on GPUs), and works directly in real space, enabling substantial wall-time reductions. The CBS-RAD method attains comparable accuracy but is constrained by its reliance on coupled-cluster energies for large systems.

FN-DMC's balance of accuracy, scalability, and parallel performance makes it a compelling alternative for predicting substituent effects on activation energies and reaction enthalpies in systems where traditional methods become computationally prohibitive. As molecular complexity increases, FN-DMC stands out as a practical and robust tool for extending the reach of high-level quantum chemical modeling.

Author contributions

Timothy B. Huber: investigation (lead); methodology (lead); formal analysis (lead); validation (lead); visualization (lead); software (lead); writing – original draft (lead); conceptualization (equal); writing – review and editing (equal). Ralph A. Wheeler: supervision (lead); conceptualization (equal); writing – review and editing (equal).

Conflicts of interest

The authors have no conflicts to disclose.

Data availability

The data supporting the findings of this study are available in the article and its supplementary information (SI). Supplementary information contains details on the statistical analyses and electronic energies computed using FN-DMC. See DOI: https://doi.org/10.1039/d5cp02300a.

In addition, input and output files for the FN-DMC calculations are openly available on Zenodo at https://doi.org/10.5281/zenodo.17059835.

Acknowledgements

The authors acknowledge the Department of Chemistry and Biochemistry at Northern Illinois University for financial support of TBH and the Center for Research Computing and Data for generous allocations of computer time on their Gaea hybrid compute cluster. The authors also thank Dr. Anouar Benali for helpful discussions and Professor Michael Papka for generously providing computational resources in the early stages of the project.

References

1. PlasticsEurope, Annual production of plastics worldwide from 1950 to 2020 (in million metric tons) Statista, 2023, https://www.statista.com/statistics/282732/global-production-of-plastics-since-1950/.
2. P. Nesvadba, Radical Polymerization in Industry, John Wiley & Sons, Ltd, 2012.
3. R. Guerrero-Santos, E. Saldívar-Guerra and J. Bonilla-Cruz, Free Radical Polymerization, John Wiley & Sons, Ltd, 2013, ch. 4, pp. 65–83.
4. N. P. Truong, G. R. Jones, K. G. E. Bradford, D. Konkolewicz and A. Anastasaki, Nat. Rev. Chem., 2021, 5, 859–869.
5. M. L. Coote, Quantum Chemical Studies of RAFT Polymerization, John Wiley & Sons, Ltd, 2021, ch. 6, pp. 139–185.
6. M. Buback, R. G. Gilbert, G. T. Russell, D. J. T. Hill, G. Moad, K. F. O'Driscoll, J. Shen and M. A. Winnik, J. Polym. Sci., Part A: Polym. Chem., 1992, 30, 851–863.
7. S. Beuermann and M. Buback, Prog. Polym. Sci., 2002, 27, 191–254.
8. S. Beuermann, M. Buback, T. P. Davis, N. García, R. G. Gilbert, R. A. Hutchinson, A. Kajiwara, M. Kamachi, I. Lacík and G. T. Russell, Macromol. Chem. Phys., 2003, 204, 1338–1350.
9. M. Buback and H. Kattner, Mol. Phys., 2021, 119, e1939452.
10. K. J. Laidler and M. C. King, J. Phys. Chem., 1983, 87, 2657–2664.
11. J. Rezác and P. Hobza, J. Chem. Theory Comput., 2013, 9, 2151–2155.
12. S. Grimme and P. R. Schreiner, Angew. Chem., 2018, 57, 4170–4176.
13. C. Riplinger, B. Sandhoefer, A. Hansen and F. Neese, J. Chem. Phys., 2013, 139, 134101.
14. C. Riplinger and F. Neese, J. Chem. Phys., 2013, 138, 034106.
15. C. Riplinger, P. Pinski, U. Becker, E. F. Valeev and F. Neese, J. Chem. Phys., 2016, 144, 024109.
16. Y. Guo, C. Riplinger, U. Becker, D. G. Liakos, Y. Minenkov, L. Cavallo and F. Neese, J. Chem. Phys., 2018, 148, 011101.
17. S. Mallick, B. Roy and P. Kumar, Comput. Theor. Chem., 2020, 112934.
18. D. G. Liakos, M. Sparta, M. K. Kesharwani, J. M. L. Martin and F. Neese, J. Chem. Theory Comput., 2015, 1525.
19. D. G. Liakos, Y. Guo and F. Neese, J. Phys. Chem. A, 2019, 90.
20. I. Sandler, J. Chen, M. Taylor, S. Sharma and J. Ho, J. Phys. Chem. A, 2021, 1553.
21. T. B. Huber and R. A. Wheeler, Phys. Chem. Chem. Phys., 2024, 26, 27536–27543.
22. G. A. Petersson, T. G. Tensfeldt and J. A. Montgomery, Jr., J. Chem. Phys., 1991, 94, 6091–6101.
23. J. W. Ochterski, G. A. Petersson and J. A. Montgomery, Jr., J. Chem. Phys., 1996, 104, 2598–2619.
24. P. M. Mayer, C. J. Parkinson, D. M. Smith and L. Radom, J. Chem. Phys., 1998, 604.
25. J. M. L. Martin and G. de Oliveira, J. Chem. Phys., 1999, 111, 1843–1856.
26. A. D. Boese, M. Oren, O. Atasoylu, J. M. L. Martin, M. Kállay and J. Gauss, J. Chem. Phys., 2004, 120, 4129–4141.
27. A. Karton, E. Rabinovich, J. M. L. Martin and B. Ruscic, J. Chem. Phys., 2006, 125, 144108.
28. J. A. Pople, M. Head-Gordon, D. J. Fox, K. Raghavachari and L. A. Curtiss, J. Chem. Phys., 1989, 90, 5622–5629.
29. L. A. Curtiss, K. Raghavachari, G. W. Trucks and J. A. Pople, J. Chem. Phys., 1991, 94, 7221–7230.
30. L. A. Curtiss, K. Raghavachari, P. C. Redfern, V. Rassolov and J. A. Pople, J. Chem. Phys., 1998, 109, 7764–7776.
31. L. A. Curtiss, P. C. Redfern and K. Raghavachari, J. Chem. Phys., 2007, 126, 084108.
32. N. J. DeYonker, T. R. Cundari and A. K. Wilson, J. Chem. Phys., 2006, 124, 114104.
33. T. G. Williams, N. J. DeYonker, B. S. Ho and A. K. Wilson, Chem. Phys. Lett., 2011, 504, 88–94.
34. A. Mahler and A. K. Wilson, J. Chem. Theory Comput., 2013, 9, 1402–1407.
35. F. Ruipérez, Int. Rev. Phys. Chem., 2019, 38, 343–403.
36. W. M. Fabian, Monatsh. Chem., 2008, 139, 309–318.
37. W. K. Trygve Helgaker and D. P. Tew, Mol. Phys., 2008, 106, 2107–2143.
38. A. Karton, S. Daon and J. M. Martin, Chem. Phys. Lett., 2011, 510, 165–178.
39. J. M. Simmie and K. P. Somers, J. Phys. Chem. A, 2015, 119, 7235–7246.
40. S. Guo, J. Liu, W. Qian, W. Zhu and C. Zhang, Energy Mater. Front., 2021, 2, 292–305.
41. N. Mardirossian and M. Head-Gordon, Mol. Phys., 2017, 115, 2315–2372.
42. C. Duan, A. Nandy, R. Meyer, N. Arunachalam and H. Kulik, Nat. Comput. Sci., 2023, 38–47.
43. M. L. Coote, Macromol. Theory Simul., 2009, 18, 388–400.
44. E. I. Izgorodina, M. L. Coote and L. Radom, J. Phys. Chem. A, 2005, 109, 7558–7566.
45. E. I. Izgorodina, D. R. B. Brittain, J. L. Hodgson, E. H. Krenske, C. Y. Lin, M. Namazian and M. L. Coote, J. Phys. Chem. A, 2007, 111, 10754–10768.
46. I. Y. Zhang, J. Wu, Y. Luo and X. Xu, J. Chem. Theory Comput., 2010, 6, 1462–1469.
47. N. Mardirossian and M. Head-Gordon, J. Chem. Phys., 2018, 148, 241736.
48. M. L. Coote, J. Phys. Chem. A, 2005, 109, 1230–1239.
49. M. L. Coote and D. J. Henry, Macromolecules, 2005, 38, 1415–1433.
50. G. Li Manni, R. K. Carlson, S. Luo, D. Ma, J. Olsen, D. G. Truhlar and L. Gagliardi, J. Chem. Theory Comput., 2014, 10, 3669–3680.
51. L. Gagliardi, D. G. Truhlar, G. Li Manni, R. K. Carlson, C. E. Hoyer and J. L. Bao, Acc. Chem. Res., 2017, 50, 66–73.
52. A. M. Sand, K. M. Kidder, D. G. Truhlar and L. Gagliardi, J. Phys. Chem. A, 2019, 123, 9809–9817.
53. E. C. Mitchell, T. R. Scott, J. J. Bao and D. G. Truhlar, J. Phys. Chem. A, 2022, 126, 8834–8843.
54. R. K. Carlson, G. Li Manni, A. L. Sonnenberger, D. G. Truhlar and L. Gagliardi, J. Chem. Theory Comput., 2015, 11, 82–90.
55. R. K. Carlson, D. G. Truhlar and L. Gagliardi, J. Chem. Theory Comput., 2015, 11, 4077–4085.
56. S. O. Odoh, G. L. Manni, R. K. Carlson, D. G. Truhlar and L. Gagliardi, Chem. Sci., 2016, 7, 2399–2413.
57. R. M. Martin, L. Reining and D. M. Ceperley, Stochastic Methods, Cambridge University Press, 2016, ch. 22–24, pp. 577–643.
58. J. Toulouse, R. Assaraf and C. J. Umrigar, Introduction to the Variational and Diffusion Monte Carlo Methods, Elsevier, 2016, ch. 15, vol. 73, pp. 285–314.
59. T. A. Anderson, M. C. Per and C. J. Umrigar, J. Chem. Phys., 2024, 160, 104110.
60. M. Dubecký, P. Jurecka, R. Derian, P. Hobza, M. Otyepka and L. Mitas, J. Chem. Theory Comput., 2013, 9, 4287–4292.
61. A. Benali, L. Shulenburger, N. A. Romero, J. Kim and O. A. von Lilienfeld, J. Chem. Theory Comput., 2014, 10, 3417–3422.
62. A. Ambrosetti, D. Alfè, R. A. J. DiStasio and A. Tkatchenko, J. Phys. Chem. Lett., 2014, 5, 849–855.
63. S. Grimme, Chem. – Eur. J., 2012, 18, 9955–9964.
64. Y. S. Al-Hamdani, P. R. Nagy, A. Zen, D. Barton, M. Kállay, J. G. Brandenburg and A. Tkatchenko, Nat. Commun., 2021, 12, 3927.
65. X. Zhou and F. Wang, J. Comput. Chem., 2017, 38, 798–806.
66. K. Krongchon, B. Busemeyer and L. K. Wagner, J. Chem. Phys., 2017, 146, 124129.
67. F. Fracchia, C. Filippi and C. Amovilli, J. Chem. Theory Comput., 2013, 9, 3453–3462.
68. A. C. Kollias, D. Domin, G. Hill, M. Frenklach, D. M. Golden and W. A. Lester Jr., Int. J. Chem. Kinet., 2005, 37, 583–592.
69. A. Hou, X. Zhou, T. Wang and F. Wang, J. Phys. Chem. A, 2018, 122, 5050–5057.
70. N. Nemec, M. D. Towler and R. J. Needs, J. Chem. Phys., 2010, 132, 034111.
71. F. R. Petruzielo, J. Toulouse and C. J. Umrigar, J. Chem. Phys., 2012, 136, 124116.
72. M. C. Per, E. K. Fletcher, E. T. Swann and D. M. Cleland, J. Comput. Chem., 2020, 41, 2378–2382.
73. X. Zhou, H. Zhao, T. Wang and F. Wang, Phys. Rev. A, 2019, 100, 062502.
74. T. B. Huber and R. A. Wheeler, J. Chem. Phys., 2024, 161, 034303.
75. H. Fischer and L. Radom, Angew. Chem., Int. Ed., 2001, 40, 1340–1371.
76. A. D. Becke, J. Chem. Phys., 1993, 98, 1372–1377.
77. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
78. S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys., 1980, 58, 1200–1211.
79. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
80. T. H. Dunning Jr., J. Chem. Phys., 1989, 90, 1007–1023.
81. J. Toulouse and C. J. Umrigar, J. Chem. Phys., 2007, 126, 084102.
82. J. Kim, A. D. Baczewski, T. D. Beaudet, A. Benali, M. C. Bennett, M. A. Berrill, N. S. Blunt, E. J. L. Borda, M. Casula, D. M. Ceperley, S. Chiesa, B. K. Clark, R. C. Clay, K. T. Delaney, M. Dewing, K. P. Esler, H. Hao, O. Heinonen, P. R. C. Kent, J. T. Krogel, I. Kylänpää, Y. W. Li, M. G. Lopez, Y. Luo, F. D. Malone, R. M. Martin, A. Mathuriya, J. McMinis, C. A. Melton, L. Mitas, M. A. Morales, E. Neuscamman, W. D. Parker, S. D. P. Flores, N. A. Romero, B. M. Rubenstein, J. A. R. Shea, H. Shin, L. Shulenburger, A. F. Tillack, J. P. Townsend, N. M. Tubman, B. V. D. Goetz, J. E. Vincent, D. C. Yang, Y. Yang, S. Zhang and L. Zhao, J. Phys.: Condens. Matter, 2018, 30, 195901.
83. Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov and G. K. L. Chan, J. Chem. Phys., 2020, 153, 024109.
84. Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters and G. K. L. Chan, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1340.
85. Q. Sun, J. Comput. Chem., 2015, 36, 1664–1671.
86. M. Burkatzki, C. Filippi and M. Dolg, J. Chem. Phys., 2007, 126, 234105.
87. M. Caffarel, T. Applencourt, E. Giner and A. Scemama, J. Chem. Phys., 2016, 144, 151103.
88. T. J. Lee and P. R. Taylor, Int. J. Quantum Chem., 1989, 36, 199–207.
89. M. Casula, S. Moroni, S. Sorella and C. Filippi, J. Chem. Phys., 2010, 132, 154113.
90. T. Zytowski and H. Fischer, J. Am. Chem. Soc., 1996, 118, 437–439.

---