A Combined Experimental and Theoretical Study of ESR Hyperfine Coupling Constants for N,N,N’,N’-Tetrasubstituted p-Phenylenediamine Radical Cations

A test set of N,N,N’,N’-tetrasubstituted p-phenylenediamines are experimentally explored using ESR (electron spin resonance) spectroscopy and analysed from a computational standpoint thereafter. This computational study aims to further aid structural characterisation by comparing experimental ESR hyperfine coupling constants (hfccs) with computed values calculated using ESR-optimised “J-style” basis sets (6-31G(d,p)-J, 6-31G(d,p)-J, 6-311++G(d,p)-J, pcJ-1, pcJ-2 and cc-pVTZ-J) and hybrid-DFT functionals (B3LYP, PBE0, TPSSh, ωB97XD) as well as MP2. PBE0/6-31g(d,p)-J with a polarised continuum solvation model (PCM) correlated best with the experiment, giving an R2 value of 0.8926. A total of 98% of couplings were deemed satisfactory, with five couplings observed as outlier results, thus degrading correlation values significantly. A higher-level electronic structure method, namely MP2, was sought to improve outlier couplings, but only a minority of couples showed improvement, whilst the remaining majority of couplings were negatively degraded.


Introduction
Electron spin resonance (ESR)-also commonly known as electron paramagnetic resonance (EPR) spectroscopy-is an important tool in the study and characterisation of organic radical species [1]. Computational investigations aiding such ESR experiments are critical in assessing especially short-lived species that are experimentally difficult to monitor.
A chemically interesting group of molecules that are an exception to this fact are persistent N,N,N',N'-tetrasubstituted p-phenylenediamine radical cationic species, the redox properties of which have been explored in detail [2,3]. Applications include materials for solar cells [4,5], artificial photosynthesis [6][7][8], and organic molecular conductors [9], to name a few. The first study of N,N,N',N'-tetrasubstituted p-phenylenediamines included ESR electron nuclear double resonance spectroscopy (ENDOR) measurements complemented by density functional theory (DFT) calculations of TMePD, TiPrPD, and BPyrB [10] (For structures and abbreviations, see Figure 1 and Table 1). At the B3LYP/6-31G(d,p) level of theory, hyperfine coupling constants (hfccs) were underestimated by 5-8%, with progress stunted since. Recently, this molecular test set was expanded, and the redox properties [2,11] have been studied in great detail. Substantial progress has been made in the measurement of formal redox potentials with computational analysis verifying experimental potentials with high confidence. Different computational methods have displayed varying degrees of success, with a multi-level G3MP2B3 method producing the highest R 2 value of 0.9624 but at the cost of a diminished slope and y-intercept. DFT/B3LYP with a large 6-311++G(d,p) basis set improved the slope and y-intercept parameters significantly, with a marginally depreciated R 2 value in comparison to G3MP2B3. TPSSh/6-311++G(2d,p) yielded the most robust results with a slightly less than perfect slope. From this, DFT appears to be the obvious choice for this articles' computational explorations. Additionally, since the initial investigation by Grampp et al. [10], ESR/NMR optimised J-augmented basis sets have become widely available [12][13][14][15][16][17][18][19], hence improvement of correlation with experimental and calculated couplings is promising.  Table 1. In this present work, we present experimental ESR parameters including 10 g tensor values and 64 hfccs for 11 N,N,N',N'-tetrasubstituted p-phenylenediamines obtained by in situ EPR/UV-vis-NIR spectroelectrochemical measurements. In addition, we report hfcc results for a series of DFT calculations aiming to determine a computational protocol suitable for reproducing experimental hfccs (similar to our previous study of N,N,N',N'-tetrasubstituted p-phenylenediamine redox potentials [2]). In principle, N,N,N',N'-tetrasubstituted p-phenylenediamines are optimal for ESR experimental investigations [20] as ESR-signals are observed only for radical species.

Experimental Results
All experimentally determined ESR hfccs and g-factors are listed in Table 2 (with DMeAzirA g-factor indeterminable). Experimental ESR results determine a H arom couplings as a positive value; yet, when compared to ENDOR investigations of similar compounds [21], a negative value is deduced. Therefore, negative experimental a H arom coupling values will be used in the linear regression studies to correctly align with computed a H arom couplings. As expected, the g-factor is invariant towards the alkyl substituents, and the value 2.0059 ± 0.0001 is obtained for all compounds. The hfccs for TMePD are a N = 0.7156 mT, a H arom = 0.2071 mT, and a H CH 3 = 0.6686 mT-differing by +0.0157 mT, +0.0074 mT, and -0.0084 mT, respectively, from reported values in similar literature [10,22,23]. Similar agreements with the literature values are observed for the hfccs of TiPrPD and BPyrB [10]. For each asymmetric compound, two nitrogen atoms are chemically different, as are the hydrogen atoms attached to the central ring, yet it is not possible to experimentally distinguish each individual atom. However, only small changes in each hfcc are observed for these atoms depending on the substituents, and their exact identification is not crucial for this article.

Computational DFT Results
In the following, the main form of statistical analyses for correlating experimental and computed results will be communicated via correlation plots and calculation of suitable quantitative measures such as R 2 values, slopes (m), and y-axis intercepts (c). The results for the individual molecules, couplings, functionals and basis sets are shown in Tables S1-S11 of the Supplementary Material.
Initially, a study using B3LYP and a relatively small, double-zeta Pople basis set, 6-31G(d,p), was conducted (see Figure 2), yielding R 2 = 0.8849, slope = 0.914, and y-intercept = 0.029 (Table 3) between the calculated and experimental couplings within the collated set of eleven target molecules. Considering that we have included a diverse set of nuclei, the correlation with experiment is relatively satisfactory.

of 17
All experimentally determined ESR hfcc's and g-factors are listed in Table 2 (with 120 DMeAzirA g-factor indeterminable). Experimental ESR results determine a H arom couplings 121 as a positive value -yet, when compared to ENDOR investigations of similar compounds 122 [40], a negative value is deduced. Therefore, negative experimental a H arom coupling values 123 will be used in the linear regression studies to correctly align with computed a H arom couplings. 124 As expected, the g-factor is invariant towards the alkyl substituents, and the value 2.0059 ± 125 0.0001 is obtained for all compounds. The hfcc's for TMePD are a N = 0.7156 mT, a H arom = 126 0.2071 mT, and a H CH3 = 0.6686 mT -differing by +0.0157 mT, +0.0074 mT, and -0.0084 mT 127 respectively, from reported values in similar literature [10,41,42]. Similar agreements with 128 literature values are observed for hfcc's of TiPrPD and BPyrB [10]. For each asymmetric 129 compound, two nitrogen atoms are chemically different, as are the hydrogen atoms attached 130 to the central ring, yet it is not possible to distinguish experimentally each individual atom. 131 However, only small changes in each hfcc are observed for these atoms depending on the 132 substituents -and their exact identification is not crucial for this article.

134
In the following the main form of statistical analyses for correlating experimental and 135 computed results will be communicated via correlation plots and calculation of suitable 136 quantitative measures such as R 2 values, slopes (m) and y-axis intercepts (c). The results for 137 the individual molecules, couplings, functionals and basis sets are shown in Tables S1 to 138 S11 of the Supplementary Material. It is quite clear that 6-31G(d,p) still reigns superior over the augmented "J-style" basis sets, retaining a higher R 2 and slope combined with an impressively lower y-intercept value. Comparing the best performing "J-style" 6-31G(d,p)-J basis set with the dominant 6-31G(d,p), one observes the largest dis-improvement in terms of slope and y-intercept as we transition to the "J-style" regime. With respect to the R 2 value, an extra outlier coupling is observed (six couplings in Figure 2    In terms of the "J-style" basis set performance, 6-31G(d,p)-J fairs best, with a degrading correlation trend observed as basis set size increases towards cc-pVTZ-J. Comparing 6-31G(d,p)-J and cc-pVTZ-J, one observes that all coupling groups become less correlated contributing to a noticeably worse R 2 value. In contrast, a better slope is observed for cc-pVTZ-J in comparison to 6-31G(d,p)-J, yet such behaviour is attributed to a larger set of under-calculated couplings, namely a N couplings, inadvertently shifting the slope towards unity. Comparing absolute calculated coupling values, 59% of nitrogen couplings become underestimated by 10% or more whilst all other coupling groups remain consistent throughout the "J-style" series, including outlier couplings. Similar trends are followed for basis set series comparisons throughout this article.
Another potential route for improvement follows changing the DFT functional and trialling each functional (TPSSh and PBE0) with a combination of each aforementioned "J-style" basis set. Comparisons between newly tested functionals and "J-style" basis sets with the energetically optimised dominant 6-31G(d,p) basis set are shown in Table 4. causal explanation of these trends eludes that the under-estimation of, namely a n couplings, 208 produces seemingly opposing trends. As shown, the results are not satisfactorily intuitive -and saying 6-31G(d,p) is a 210 "better" basis set than cc-pVTZ-J is incorrect. One potential source of error within this 211 works' computational model lies with the implementation of an implicit solvation model 212 specifically the polarised contiuum model (PCM). It is hard to judge how much this 213 may effect the accuracy of the calculation but for some systems this may be considerable. 214 Instead one could implement an explicit solvation model to better replicate an experimental 215 solvation envionment. As shown in Table 8, the improvement is considerable transitioning 216 from a system absent of solvation to a system with implicit solvation, hence this more 217 accurate model may bring potential improvements. In reality, due to the high computational 218 cost associated with an explicit solvation model paired with this studies relatively large 219 molecules, an explicit route may not be advised.   With reference to Table 4, the best performing basis set and functional combination is PBE0/6-31G(d,p)-J/PCM with an R 2 = 0.8926, m = 0.873 and c = 0.039. Arguably, a higher R 2 value of 0.8957 (vs. 0.8926) is recorded for TPSSh/6-31G(d,p)-J/PCM, but the strength of correlation (i.e., slope) diminished (m = 0.845 vs. 0.873) considerably in comparison to PBE0/6-31G(d,p)-J/PCM. The relative error (i.e., intercept) for each model is equal at c = 0.039. Therefore, PBE0/6-31G(d,p)-J/PCM is deemed a superior model due to a noticeably better strength of correlation. The relative lack of description within the 6-31G(d,p)-J basis leads one to conclude that error cancellation effects may be of total dominance in this study. The inherent issue may arise from the error of the approximated E XC [ρ] functional cancelling with the large, direct basis set error of the smallest 6-31G(d,p)-J basis set.

Discussion
Literary articles [10,24] disseminating the calculation of hfccs typically employ small basis sets in hybrid-DFT-ESR calculations, concluding that an error cancellation approach may be counter-intuitively optimal in this instance. With reference to TPSSh and B3LYP, the gap between the worst and best described basis sets, 6-31G(d,p)-J and cc-pVTZ-J, respectively, is more evident, with PBE0 actually maintaining a more negligible gap between 6-31G(d,p)-J and cc-pVTZ-J. In retrospect, the PBE0 functional best abides with a logical a priori prediction of the largest cc-pVTZ-J basis set best aligning calculated results with experimental coupling constants. From this, one may conclude that the PBE0 functional may suffer least from functional error by retaining the most satisfactory cc-pVTZ-J results.
A marginal improvement of ≈1% is seen in Figure 3 using PBE0/6-31G(d,p)-J/PCM model for calculating a H CH 2α and a H CH 2 β couplings in comparison to B3LYP/6-31G(d,p), hence promoting a greater R 2 value. However, outlier couplings degraded further, giving a marginally poorer correlation slope and y-intercept as shown in Table 4. It is difficult to conclude which results are better between each model, but based on the hierarchical importance of the R 2 value, the PBE0/6-31G(d,p)-J/PCM model is subjectively better. As seen in Table 3, two opposing trends also appear here in Table 4 as we correlate the change in basis set size with the linearly regressed parameters; R 2 , slope and y-intercept across the three functionals. Trending the change in R 2 and y-intercept, one observes improved values for each parameter as basis set size shrinks. Conversely, when examining the change in slope, a value closer to unity is calculated as the basis set size increases. The fact that underestimation of, namely a n couplings, produces seemingly opposing trends eludes a causal explanation of these trends.
As shown, the results are not satisfactorily intuitive-and saying 6-31G(d,p) is a "better" basis set than cc-pVTZ-J is incorrect. One potential source of error within this works' computational model lies with the implementation of an implicit solvation model, specifically the polarised continuum model (PCM). It is hard to judge how much this may effect the accuracy of the calculation but for some systems this may be considerable. Instead, one could implement an explicit solvation model to better replicate an experimental solvation environment. As shown in Table 5, the improvement is considerable transitioning from a system absent of solvation to a system with implicit solvation, hence this more accurate model may bring potential improvements. In reality, due to the high computational cost associated with an explicit solvation model paired with this studies relatively large molecules, an explicit route may not be advised.

Basis Set Sensitivity and Robustness
It is also of interest to examine correlation parameter's sensitivity and robustness to increasing basis set size for each functional. Firstly, regarding B3LYP, R 2 and slope change by less than ±2.1%, whilst PBE0 is best, at less than ±1.5%. Meanwhile, TPSSh functional robustness is poorest with changes of less than ±3.1%, concluding that TPSSh sensitivity to change in basis set size is greatest. All functionals show quite poor robustness in terms of the y-intercept parameter; however, these values are truly satisfactory as they near zero, hence arguments of robustness may be deemed redundant.
Examining cc-pVTZ-J for all functionals, the trend in R 2 and y-intercept (except PBE0) improve when transitioning from pcJ-2, suggesting that couplings begin to recover and correlate marginally better as shown in Figure 4. Shown also is the depreciating trend in slope (except B3LYP) as we transition from pcJ-2 to cc-pVTZ-J. From this, one concludes that overall coupling correlation recovers, yet calculated couplings remain overestimated when the largest cc-pVTZ-J basis set size is realised. the y-intercept parameter -however, these value are truly satisfactory as they near zero, 227 hence arguments of robustness may be deemed redundant.
228 Figure 4. Trend in R 2 , slope (m) and y-intercept (c) parameters transitioning along each basis set for B3LYP, TPSSh and PBE0 functionals Examining cc-pVTZ-J for all functionals, the trend in R 2 and y-intercept (except PBE0) 229 improve when transitioning from pcJ-2 suggesting that couplings begin to recover and 230 correlate marginally better as shown in Figure 4. Shown also is the depreciating trend in 231 slope (except B3LYP) as we transition from pcJ-2 to cc-pVTZ-J. From this, one concludes 232 that overall coupling correlation recovers, yet calculated couplings remain overestimated 233 when the largest cc-pVTZ-J basis set size is realised. Although, common to all "J-style" basis set models tested are the constant set of 236 detrimental outliers degrading calculated coupling quality significantly -hence an analytic 237 effort is applied hereafter on diagnosing the origin of such erroneous behaviour.

238
As shown in Figure 3, problematic outliers remain drastically depreciating the quality 239 of the R 2 value (0.873), slope (0.8926) and y-intercept (0.039) within the best performing 240 "J-style" system PBE0/6-31G(d,p)-J. Regardless of basis set/functional deployed, these 241 outliers still stem from two molecules exclusively; DMeAzirA and DMeAzetA (shown 242 in Figure 5) for all "J-style" basis set calculations. From initial observations, this is an 243 Figure 4. Trend in R 2 , slope (m) and y-intercept (c) parameters transitioning along each basis set for B3LYP, TPSSh and PBE0 functionals.

Examining Outliers
Common to all "J-style" basis set models tested is the constant set of detrimental outliers degrading calculated coupling quality significantly; hence, an analytic effort is applied hereafter for diagnosing the origin of such erroneous behaviour.
As shown in Figure 3, problematic outliers remain, drastically depreciating the quality of the R 2 value (0.873), slope (0.8926), and y-intercept (0.039) within the best-performing "J-style" system PBE0/6-31G(d,p)-J. Regardless of the basis set/functional deployed, these outliers still stem from two molecules exclusively for all "J-style" basis set calculations: DMeAzirA and DMeAzetA (shown in Figure 5). From initial observations, this is an issue pertaining to strained fragments within each molecular structure that are not treated correctly by any hybrid-DFT functional and/or basis set tested.
From a structural standpoint, each molecule shares a highly strained N-ring that projects significant problems onto the coupling calculation for each group present within the ring-with couplings being either drastically over/under-estimated. Four out of five outliers are present within each strained fragment, with the fifth outlier a H CH 3 attached to the non-cyclic-nitrogen in DMeAzirA. On average, couplings appeared less correlated with experimental values in DMeAzirA compared to DMeAzetA, concluding that increased straining effects have a global impact on couplings within the molecule. In comparison, couplings within the slightly less strained DMeAzetA system are better correlated with experimental results, yet still exhibit major outlier behaviour within the collated set of calculated couplings. Other couplings that are deemed "non-outliers" in the collated set reflect that they are independent of minor perturbations of other substituents in the instances that these perturbations do not originate from a strained ring fragment.
No other N-strained ring is present in this subset of N,N,N',N'-tetrasubstituted pphenylenediamines, concluding that this is an issue of ineptness within the computational model's ability to treat couplings associated with highly strained N-ring systems. From molecular structure theory, the interior angle between atoms in three/four member rings is much less than ideal (compared to similar five or six member rings), forcing an increased overlap of electron density. In strained systems, some "wiggle" room is allowed (maybe forcing a group out of the plane) to alleviate this overlap, so a different conformation may be optimal for the DFT functional and/or basis set used.
Referencing Table 6, a clear depreciating trend in correlation parameters is shown as one transitions from three-member towards six-member ring systems when assessing individual molecular performance with PBE0/6-31G(d,p)-J. DMeAzirA exhibits the largest ring strain within the three-member ring, imposing a poor coupling correlation, whereas the four-member DMeAzetA exhibits considerably better results, albeit still performing poorly within the set of the eleven molecules tested. In contrast, BPipB showcases the best correlation, with a marginal improvement observed transitioning from a five-member ring to a six-member ring. Such results reiterates the detrimental impact that ring strain imposes on smaller (three-four) ring systems. Of interest also was examining the change in correlation values upon removal of both problematically strained systems, shown in Table 7. As expected, all regression parameters improved significantly for all functional and basis set combinations, with similar trends to those observed in Table 4. Most evident is the superior correlation trend attained by the smallest basis set-6-31G(d,p)-J, largely due to suspected basis set and functional error cancellation. As the basis set size increases towards cc-pVTZ-J, experimental and calculated values deviate due to calculated couplings becoming under-calculated, generally, with PBE0 functional retaining greatest robustness towards basis set size increase.
These strained fragments were also deemed problematic in the calculation of oxidation potentials in a related article [2], yet it was observed that a higher-level electronic structure model treated such outliers appropriately, thus improving results significantly. In this article, a similar approach is trialled to improve DFT results by using a slightly more sophisticated model such as MP2. Due to the inordinate scaling of CCSD or CCSD(T), MP2 with N 5 (where N is one electron basis function) scaling possesses a satisfactory middle ground between computation time and accuracy.

Computational MP2 Analysis
Taking inspiration from the success of the double-zeta Pople basis sets in the previous section, it was decided to re-optimise each structure using MP2/6-31G(d,p), utilising a J-augmented 6-31G(d,p)-J basis set for each MP2 single-point hfcc calculation. Three molecules were chosen for sampling the quality of MP2, namely the two outlier molecules (DMeAzetA and DMeAzirA) and one molecule that exhibited an excellent correlation with experimental results (DMeDiPrPD). A 1:1 comparison of each coupling with the blue ideal fitting line (y = x) is shown in Figure 6 to observe the improvement, if any, that is gained transitioning from PBE0 to MP2.
Comparing changes transitioning from PBE0 to MP2 for each molecule in Figure 6, it is apparent that the discrepancies between experimental and calculated values increase significantly as we compare the well-behaved DMeDiPrPD system with the ill-behaved DMeAzetA and DMeAzirA systems. This behaviour is quite possibly a direct effect of the strained rings-most evidently seen in DMeAzirA and reducing in effect as we transition to DMeAzetA-with DMeDiPrPD behaving appropriately due to an absence of such straining effects. Such large changes suggest that strained compounds are more sensitive to how electron correlation is treated compared to non-strained compounds, as one might have expected. However, the large changes attained from going to MP2 do not improve on the quality of correlation as was the case in the previous study of redox potentials related to such N,N,N',N'-tetrasubstituted p-phenylenediamines [2].
Regarding individual coupling group comparisons, MP2 always gives poorer a H aromatic couplings in comparison to PBE0, with a similar trend observed for a N couplings (except for one a N coupling in DMeAzetA). On average, a H CH 2 α couplings are split 50:50 between each model, with a H CH 2 β and a H CH 3 couplings actually improving from PBE0 to MP2. One a H CH coupling was present in favour of PBE0. From this, one can conclude that particular coupling groups are degraded when transitioning from PBE0 to MP2, but others may or may not improve. Comparing changes transitioning from PBE0 to MP2 for each molecule in Figure 300 6, it is apparent that discrepancies between experimental and calculated values increase 301 significantly as we compare the well-behaved DMeDiPrPD system with the ill-behaved 302 DMeAzetA and DMeAzirA systems. This behaviour is quite possibly a direct effect of the 303 strained rings -most evidently seen in DMeAzirA and reducing in effect as we transition to 304 DMeAzetA -with DMeDiPrPD behaving appropriately due to an absence of such straining 305 effects. Such large changes suggest that strained compounds are more sensitive to how 306 electron correlation is treated compared to non-strained compounds, as one might have 307 expected. However, the large changes attained from going to MP2 do not improve on the 308 quality of correlation as was the case in the previous study of redox potentials related to 309 such N,N,N',N'-tetrasubstituted p-phenylenediamines [2].

310
Regarding individual coupling group comparisons, MP2 always gives poorer a H aromatic 311 couplings in comparison to PBE0, with a similar trend observed for a N couplings (except 312 for one a N coupling in DMeAzetA). On average, a H CH 2 α couplings are split 50:50 between 313 each model -with a H CH 2 β and a H CH 3 couplings actually improving from PBE0 to MP2. One 314 a H CH coupling was present in favour of PBE0. From this, one can conclude that particular 315 coupling groups are degraded when transitioning from PBE0 to MP2, but others may, or 316 may not improve. The question arises now whether the changes going from PBE0 to MP2 calculations are due to changes in geometry upon performing an MP2 geometry optimization or else due to MP2 calculation of spin density at the nuclei. Hence, it is of interest to compare the hfcc results of a single-point PBE0/6-31G(d,p)-J and MP2/6-31G(d,p)-J calculation with both retaining MP2/6-31G(d,p) optimised geometries for each outlier molecule, DMeAzetA and DMeAzirA.
As shown in Figure 7, one observes that the PBE0 hyperfine couplings remain stagnant compared to Figure 6 differences between PBE0 and MP2 calculated values originating solely from the differences in the spin-density in the hfcc MP2 single-point calculation. Strained effects are still in total effect in both MP2 and PBE0 optimised systems for DMeAzetA and DMeAzirA, with the geometry change between PBE0 and MP2 contributing little to no effect towards improved hfccs (when comparing each molecules performance in Figures 6 and 7). Therefore, the following regression analysis will be based with Figure 6 in mind. The usual linear regression model is used in Table 8 to compare R 2 , slope and yintercept values retained by MP2 and PBE0. For ease of comparison, both problematic systems and a model DMeDiPrPD system (that achieved satisfactory results) are used to statistically determine if any improvement is gained transitioning from the best performing DFT functional PBE0 to a higher-level method such as MP2. Table 8. Linear regression analysis of MP2, PBE0, and ωB97XD for both problematic systems (DMeAzeta and DMeAzirA) with reference to a satisfactory model system (DMeDiPrPD). Evidently, PBE0 remains superior to MP2 for DMeAzetA and DMeDiPrPD as the R 2 , slope, and y-intercept retain better values. With respect to DMeAzirA, one actually sees an improved R 2 value, whereas a quite drastic degradation in slope and y-intercept is concurrently observed when transitioning from PBE0 to MP2. These opposing correlation trends for DMeAzirA are explained simply by observing the holistic data point distribution with respect to the fitted line in each separate plot of MP2 and PBE0 in Figure 6. A closer to linear relationship is seen with a superior R 2 for MP2. Yet, factoring a depreciated slope and y-intercept, one concludes that there is a large error from the true fitted line. Moreover, this error originates from the overestimation of calculated couplings at the MP2 level of theory, which is clearly seen upon comparison with corresponding absolute experimental values.
A reference calculation was also performed ( Table 8) using ωB97XD optimised using MP2/6-31G(d,p) and a single-point calculation ran with a 6-31G(d,p)-J basis set and PCM solvation model. Surprisingly, DMeAzirA improved significantly over MP2 and PBE0 showing signs for improved correlation upon realisation of the ωB97XD functional. Meanwhile, DMeAzetA and DMeDiPrPD correlation values remain stationary in comparison to MP2 and PBE0. In order to compare this result when functionals previously tested, it was necessary to perform a similar single-point calculation with an ωB97XD/ 6-31G(d,p) optimised geometry for each of the three test molecules. It is immediately apparent that correlation parameter performance degrades significantly for DMeAzirA in comparison to the MP2/6-31G(d,p) optimised equivalent calculation. Thus, this result falls in line with similar calculations using B3LYP, TPSSh and PBE0 functionals.
As a means to further consolidate discussion of the MP2 results, it was of interest to run a simple single-point calculation using the CFOUR [25] quantum chemistry program. An MP2 calculation was run with an MP2/6-31G(d,p) optimised geometry without solvation parameters and compared to an equivalent calculation run in the Gaussian program. CFOUR by implementation only prints the spin density of each nuclei in the output, hence an extra conversion is needed to calculate the isotropic hfcc a iso of each nuclei N: where S z is the total spin expectation value, g e and g n are the electron and nuclear gfactors, and β e and β N are the electron and nuclear magnetrons. ρ( R N ) is the nuclear spin density. As shown in Table 5, the results are significantly similar to conclude the MP2 results presented in Table 8 are consistent. Shown also are the DMeAzetA experimental results, which are in much better agreement with MP2 Gaussian results in a solvated environment, further exemplifying the critical improving effect a solvation model provides. Therefore, to calculate N,N,N',N'-tetrasubstituted p-phenylenediamines hfccs, using MP2 would be irrational as the benefit is only marginal for a subset of couplings taking into consideration also the extra cost with respect to hybrid DFT methods. The final synopsis would be to therefore continue with PBE0/6-31G(d,p)-J/PCM for calculating hfccs of N,N,N',N'-tetrasubstituted p-phenylenediamines, but a case can be made for using ωB97XD with an MP2-optimised geometry, especially for improving couplings in the strained DMeAzirA molecule.

Materials and Methods
In the experiment outlined, ESR and UV-vis-NIR measurements are recorded simultaneously, hence both will be referenced within the experimental procedure. Disseminated results from this experiment concentrate solely on ESR measurements and computed results; hence, UV-vis-NIR results are not discussed. Detailed syntheses of each compound have been previously described in detail [2]. Furthermore, compounds were kept under argon during travel to Bratislava, Slovakia and until used for the ESR experiments. The purity of the compounds was then checked by cyclic voltammetry in Bratislava, Slovakia.

Experimental Procedure
Samples for in situ ESR/UV-vis-NIR spectroelectrochemical experiments were prepared by dissolving a weighted amount of electroactive solute and supporting electrolyte (Bu 4 NPF 6 ) in argon degassed acetonitrile. The concentrations of the solute and the electrolyte were 1 mM and 0.2 M, respectively. All cyclic voltammograms were recorded under argon atmosphere in the potential range of the solute, and the potentiostat triggered the ESR and UV-vis-NIR spectrometers in order to synchronize measurements. The spectroelectrochemical cell (0.1 mm path length) was specially designed for the experiment [26] and suitable for an optical transmission ESR resonator (ER 4104 OR-C 9609) of an EMX ESR spectrometer. The working electrode was a laminated Pt mesh with a small hole in the foil coincident with the light beam, which limited the active surface area of the electrode. A Pt wire auxiliary (counter) electrode and a Ag wire pseudoreference electrode were used. The optical ESR resonator cavity was connected to the diode-array UV-vis-NIR spectrometer Avantes Avaspec (Avantes, Apeldoorn, The Netherlands) by optical fibres. A deuterium-halogen lamp DH 2000 (Sentronic, Schleswig-Holstein, Germany) was used as a light source. UV-vis-NIR spectra were processed by the proprietary AvaSoft 7.7 software package, and the experimental ESR spectra were analysed and simulated using the Bruker softwareWinEPR and SimFonia, respectively. The g-values of the radical cations were determined by using MgO (Mn 2+ ) as an internal standard.

Computational Procedure
Standard one-electron energetically optimised Gaussian basis sets are, in general, not well suited for calculation of hfccs. Conceptually, ESR-optimised basis sets aim to better quantify the true electron density near the nucleus, better approximating the isotropic Fermi-contact term as a result. The Fermi-contact mechanism acts in a coupling regime between electronic spin and the magnetic dipole of the nucleus-with the observed splitting distance between peaks in spectra defined as the magnitude of the hfcc. An anisotropic spin-dipolar interaction also constitutes part of this interaction but is averaged to zero in a solvated environment.
Achieving a higher resolution of the electronic density around the nucleus through the use of very-high-exponent s-functions within each basis set is quite straight forward. Yet, true replication of experimental couplings still eludes researchers as other factors such as electron correlation, basis set size, or even relativistic effects always complicate matters. Furthermore, computing discrete solvation effects [42] is neither cheap or trivial, with vibrational averaging [43] also non-negligible in some cases. The latter was not analysed as a method for how to replicate solvation effects in a vibrationally averaged regime has yet to be disseminated. Instead, a polarised continuum solvation model (PCM) provided a satisfactory middle point. Following this, the mantra of this study was not to replicate experimental results with total accuracy, but to provide a clear rationalization of hfcc trends and link outlier results to possible structural anomalies.
Computationally, each target structure underwent a full geometrical optimisation with the same computational model as the calculation of the hfccs with the exception of using standard energetically optimised basis sets instead, i.e., the "non-J augmented" forms. Afterwards, hfccs were calculated within a single-point calculation using the "J-augmented" version of each basis set. In every calculation, the polarised continuum (PCM) solvation model (Gaussian09 keyword; IEFPCM) was deployed replicating fundamental solvation in an approximated framework. All calculations were carried out with the GAUSSIAN09 [44] program. All hfccs are given in milliteslas, mT (1 mT = 10 Gauss).

Conclusions
In essence, we have presented experimental values of the ESR hyperfine coupling constants of 11 structurally related N,N,N',N'-tetrasubstituted p-phenylenediamines. As a reference, we have tested the performance of different computational protocols using a large set of 64 hyperfine couplings consisting of 17 nitrogen and 47 hydrogen nuclei. We concentrated mostly on DFT utilising three exchange correlation functionals-B3LYP, PBE0, and TPSSh-and also presented MP2 calculations for selected systems. Furthermore, we have investigated the performance of "J-style" basis sets which are specially optimised for the calculation of coupling constants.
Overall, calculated values correlated well with experimental values in the scenario that outlier DMeAzetA and DMeAzirA hfccs were excluded from the regression analysis. The usage of "J-style" ESR-optimised basis sets also retained varying degrees of success, with the general trend of increasing basis set size counter-intuitively invoking a degradation of hfcc quality as we moved from a small 6-31G(d,p)-J basis set towards the largest cc-pVTZ-J set. Amongst the set of hybrid DFT functionals tested, PBE0 performed marginally best, hence uniting PBE0 and 6-31G(d,p)-J as the best performing hybrid DFT functional and basis set, respectively, within a PCM solvated environment. Moreover, a transition to MP2 as the electronic structure method yielded a large change in results, especially for both strained systems DMeAzetA and DMeAzirA, yet in a direction that provided no discernible improvement. Hence, PBE0/6-31G(d,p)-J is the choice of electronic structure method and basis set, respectively, within a PCM solvated environment for this set of 11 N,N,N',N'-tetrasubstituted p-phenylenediamines.

Data Availability Statement:
The data that support the findings of this study are available within the article and the Supporting Information. Further data are available from the corresponding author upon reasonable request.