A Computational Journey across Nitroxide Radicals: From Structure to Spectroscopic Properties and Beyond

Nitroxide radicals are characterized by a long-lived open-shell electronic ground state and are strongly sensitive to the chemical environment, thus representing ideal spin probes and spin labels for paramagnetic biomolecules and materials. However, the interpretation of spectroscopic parameters in structural and dynamic terms requires the aid of accurate quantum chemical computations. In this paper we validate a computational model rooted into double-hybrid functionals and second order vibrational perturbation theory. Then, we provide reference quantum chemical results for the structures, vibrational frequencies and other spectroscopic features of a large panel of nitroxides of current biological and/or technological interest.


Introduction
Nitroxides are one of the few families of stable organic free radicals and this feature, together with the remarkable sensitivity of their structure and spectroscopic properties to environmental effects has stimulated their widespread use as spin labels and spin probes in both biological and material chemistry [1][2][3]. The interest for this class of radicals has seen very recently a remarkable increase especially in connection with the analysis of dynamical and environmental effects by state-of-the-art computational approaches and with protein crystallography refinement [4][5][6]. Moreover, while the EPR (electronic paramagnetic resonance) spectra of a huge number of molecular systems including the NO moiety have been recorded and interpreted by means of quantum chemical (QC) computations [7][8][9][10], the situation is different concerning accurate molecular structures and, especially, vibrational spectra.
As a matter of fact, the EPR studies of spin labelled proteins can be connected with crystallographic experiments allowing the determination of a static atomistic description of the whole molecular model through structure refinement procedures [6]. The refinement requires additional a priori information, supplied in the form of chemical restraints, to compensate for the lack of high-resolution data and/or other experimental issues. Of course, the employed restraints should be at least as accurate as the sought accuracy of the refinement and can be often derived from accurate experimental data available for suitable fragments of the investigated molecular system. However, this approach is seldom exploitable for open-shell moieties, whose accurate structures are usually unknown. Under such circumstances QC computations can be effectively employed to obtain the missing information, provided that the accuracy of the selected computational model is sufficient [11]. In recent years, methods rooted in the density functional theory (DFT) have emerged as the methods of choice for medium-to large-size molecular systems not amenable to the most accurate (but also prohibitively expensive) wave-function methods. In this connection, several studies have shown that double-hybrid functionals in conjunction with partially augmented triple-zeta basis sets are particularly reliable for geometrical structures, vibrational frequencies and several other spectroscopic properties [7,12,13]. On these grounds, we have performed a comprehensive study of several prototypical nitroxide radicals using as benchmarks either available experimental data or state-ofthe-art QC results. After validating the selected computational approach we provide reference data for a large panel of nitroxides of current interest for biological and/or technological applications.

Computational Details
Unless explicitly stated, calculations were performed with the GAUSSIAN16 suite of programs [14]. The underlying electronic-structure model is rooted in the density functional theory (DFT). The hybrid B3LYP [15,16] functional in conjunction with the SNSD basis set [17] (hereafter B3) and the double hybrid B2PLYP [18] functional in conjunction with the maug-cc-pVTZ basis set [19] (where d functions on hydrogen atoms have been removed) (hereafter B2) were employed. Even if not explicitly indicated, Grimme's D3 empirical dispersion corrections [20] with the Becke-Johnson damping [21] (D3BJ) were always used. In order to mimic experimental conditions, when required, bulk solvent effects were taken into account by means of the polarizable continuum model [22] within its integral equation formalism (IEF-PCM) [23]. Geometry optimizations were carried out employing very tight convergence criteria, and minima were confirmed by Hessian evaluations. Harmonic force fields were obtained using analytic energy derivatives, while higher-order derivatives were computed through numerical differentiation using a step of 0.01 √ amu a 0 for the displacements along the mass-weighted normal coordinates. VTP2 calculations relied on the GVPT2 approach, taking possible resonances into the proper account [11,17].
The potential energy surface (PES) ruling the large amplitude inversion motion around the nitrogen atom of the nitroxide moiety was sampled by relaxed scans and then used to solve numerically a one-dimensional Schrodinger equation by means of a quasi-variational approach rooted in the discrete variable method [24][25][26][27][28].
Isotropic hyperfine coupling constants (hcc) of free radicals observed in ESR spectroscopy are determined by the electron spin density at or near the position of any magnetic nucleus. In particular, the 14 N hcc (the only one considered in the following and referred to as a N ) is obtained by multiplying the corresponding spin density by 323.13 for obtaining data in MHz (115.3 for results in Gauss). Since hcc's are particularly sensitive to both very tight and diffuse s functions, the purposely tailored EPR(III) basis set has been used for these computations [8].
Finally, in order to check the accuracy of the geometries obtained at the B2 level, the structure of the prototypical dimethylnitroxide radical (a) was optimized also at the CCSD(T)-F12/cc-pVDZ-F12 level Supplementary Materials, which is known to provide remarkably accurate results [29]. These computations were performed by the MOLPRO package [30]. In the same vein reference a N values were computed at the CCSD(T) level employing the CFOUR program [31].

Results and Discussion
The structure and labeling of the nitroxides investigated in the present paper are sketched in Figure 1. The key geometrical parameters ruling the nitroxide properties are the NO bond length (mainly affecting the magnetic g tensor and the UV-vis absorption) and the pyramidality of the nitrogen atom (mainly affecting the hyperfine tensor). Noted is that both quantities (which are not fully unrelated) are sensitive to the overall molecular backbone.
On these grounds, nitroxides can be classified in three principal groups: acyclic molecules (a-f); cyclic species containing six-(g-i) or five-(j-m) member rings, with the latter family being further split in two by the lack (j,k) or presence (l,m) of unsaturated bonds within the cycle. More complex systems containing aromatic moieties (n-p) have been also taken into account.

Template Molecule
It is well documented that B2 equilibrium geometries are remarkably accurate [32]; nonetheless higher level calculations performed on small model systems can provide information to correct the structures of larger molecules obtained at lower computational levels following the template molecule and/or linear regression approaches (referred to in the following as TMA and LRA, respectively) [13].
In the present connection, dimethyl nitroxide (a) is a suitable reference system. The geometries of the minimum and the planar transition state were optimized using the CCSD(T)-F12 ansatz [33] in conjunction with the cc-pVDZ-F12 basis set [34] (hereafter DZF12). As already mentioned, this level of theory leads to results comparable with those delivered by conventional CCSD(T) computations employing much larger basis sets [29].
In Table 1 the geometrical parameters of the a molecule optimized at different levels of theory are reported. Confirming previous observations for the NO bond lengths in closed-shell systems [35], despite their limited computational cost, B2 and B3 results are remarkably accurate, with the largest deviation with respect to the CCSD(T)-F12 reference being 0.003 Å. Also the CNC valence angle is well reproduced with errors below 0.3 degrees with respect to the CC reference for both B3 and B2 methods. However, the B2 method clearly outperforms the B3 counterpart for NC bond lengths, reducing the error from 0.005 to 0.001 Å. These results suggest that the geometries of larger nitroxides optimized at the B2 level can be considered sufficiently reliable, thus not requiring any ad hoc correction. The description of the out-of-plane angle is more demanding concerning both the equilibrium value (overestimated by 3 and 2.5 degrees at the B3 and B2 level, respectively) and the height of the energy barrier to planarity (which increases from 287 to 334 and 404 cm −1 when going from B3 to B2 and CCSD(T)-F12 computations). It is well known that inversion motions ruled by low barriers to planarity are generally ill described by harmonic models [28], so that the effective structure does not coincide with the energy minimum even at low temperatures. Thus, this large amplitude motion (LAM) was investigated through a relaxed scan along the NO out-of-plane coordinate and it was used to solve numerically the vibrational Schrödiger equation with the help of a quasi-variational approach employing the discrete variable representation (DVR) [24][25][26][27][28].
In Figure 2 the PES along the LAM is reported together with the associated vibrational levels obtained with the DVR method and the corresponding vibrational wave functions of those states relevant at 300 K (population above 0.5%). By effect of the LAM, the vibrationally averaged value of the out of plane angle decreases by about 2°with respect to the equilibrium value. DFT methods are known to underestimate the absolute values of a N in nitroxide systems [4]. On these grounds, the performance of the double hybrid B2PLYP functional in conjunction with the purposely tailored EPR(III) basis set was evaluated for the Me 2 NO · model system taking the results obtained at the CCSD(T) level (hereafter CC) as references. The results collected in Table 2 show that the B3 results consistently underestimate the CC reference values (by about 5 MHz) and that the B2 method represents a remarkable improvement leading to quantitative agreement with CC for the pyramidal equilibrium structure and a slight underestimation (about 1 MHz) for the planar transition state. This trend is not surprising since for planar structures a N is determined only by spin-polarization (which is particularly hard to be described due to spin contamination and role of exact exchange). As a matter of fact, increasing the contribution of Hartree-Fock exchange (cfr. B3 and BHLYP results in Table 2) leads to better results, which are further improved by inclusion of a fraction of second order many-body perturbation (B2 results). Furthermore, Figure 3 shows that the CC, B2, and B3 results as a function of θ (the angle formed by the NO bond with the CNC plane) are quantitatively fitted by the function with a = −71.65, −78.78, −78.53 and b = 102.35, 108.53, 103.47 MHz for CC, B2, and B3 methods, respectively. This is exactly the expected behaviour since θ drives the increase of the contribution of nitrogen s orbitals to the π SOMO (singly occupied molecular orbital) of nitroxides.   These results can be employed to correct the a N for larger nitroxides by a composite scheme in which the difference between CC and B2 a N is estimated by the respective values issuing from Equation (1) at a θ angle equal to that of the target molecule.
Finally, vibrational averaging along the inversion LAM decreases the value of a N at 300 K by 1.3, 1.6 and 2.5 MHz at the B3, B2 and CC level of theory, respectively. The corresponding contribution obtained from the perturbative VPT2 approach is 1.8 MHz, showing a remarkable agreement with the variational counterpart.

Geometries
The geometrical parameters describing the C 1 C 2 NO moiety of all the nitroxides shown in Figure 1 computed at the B2 level are collected in Table 3.
In the case of the heterogeneous class of acyclic molecules, the geometrical parameters strongly depend on the nature of the substituent. It is quite apparent that bulky substituents decrease the pyramidality of the nitrogen atom leading to nearly-planar structures, which become exactly planar in the presence of aryl substituents due to more effective delocalization of π electrons.
As already pointed out in previous studies [9,36,37], a nearly planar geometry of the NO group is generally observed when it is included in a five-member ring, which becomes pyramidal when included in a six-member ring. All the six-member ring molecules share the piperdine-N-oxyl scaffold, which usually adopts a chair-like structure with the NO group occupying an equatorial position. In the saturated five-member ring systems, only small deviations from the planarity of the NO nitrogen are observed, with the ring displaying a twisted structure with the atoms opposed to the nitrogen lying outside the mean plane of the system. Resembling the case of the aryl substituents, in the acyclic system the presence of an aromatic ring forces the planarity of the nitroxides. Nonetheless, when the substituents on the double bond are not strongly conjugated, as in the case of l, small distortions of the rings are observed. In terms of Cremer and Pople coordinates [38], the total puckering amplitude is small (below 0.1 Å) if compared to the values generally observed in other systems containing five member rings (usually between 0.2 and 0.35 Å) [35]. Although this leads to an out-of-plane angle of about 10°in the most stable structure, a scan of the potential energy surface along the out-of-plane coordinate revealed that the planar structure is easily accessible even at very low temperatures. Noted is that that two equivalent minima with opposite signs of the NO out-of-plane angle exist and that their interconversion requires also the mirror conformation of the side chain.
The NO bond length shows little variability for all cyclic nitroxides, with its value being shorter by about 0.01 Å for five-member rings. With the exception of b, which bears two strongly withdrawing CF 3 groups, acyclic nitroxides have longer NO bond lengths ranging from 1.280 to 1.286 Å.

NO Stretching Frequencies
Despite the ubiquitous presence of the NO group in radical traps and EPR spin probes, little attention has been paid to its vibrational characterization. At variance with other functional groups like (e.g., CO or CN), the infrared (IR) intensity of the NO stretching is generally low and its identification is not always straightforward, thus leading to inconsistent assignments [39]. Recently, Rintoul and coworker [39] reviewed the available experimental vibrational spectra of nitroxides and performed a comprehensive comparison with the results issuing from DFT calculations. A scaling factor of 0.976 was applied to DFT results in order to account for systematic errors related to the limits of the computational model and/or to the underlying harmonic approximation. Aided by the computational results, the authors found several incorrect band assignments and proposed the range between 1340 and 1450 cm −1 for the ν NO vibration.
Within the set of nitroxide radicals collected in ref. [39] for which experimental results are available, we selected five acyclic nitroxides, together with seven cyclic radicals (four five-member and three six-emember rings), whose experimental data are sufficiently reliable. When the experimental spectra were recorded in solution, bulk solvent effects were taken into account by mean of the polarizable continuum model [23]; nonetheless, most of the experimental values were recorded in KBr or Nujol, whose perturbing effect was neglected in the computations. The NO stretching frequencies computed at B3 and B2 levels for all the nitroxides sketched in Figure 1 are compared to the available experimental data in Table 4.
It is well known that harmonic calculations overestimate the vibrational frequencies, so that empirical scaling factors are usually employed to improve the agreement with experiment. When the dimension of the system does not allow full anharmonic calculations by high-level methods, an effective alternative not involving any empirical factor is offered by the so-called hybrid approaches in which harmonic frequencies at a high level of theory (here B2) are coupled to anharmonic corrections evaluated at a lower level (here B3) leading to ν H = ω B2 + (ν B3 − ω B3 ). The reliability of the issuing B2/B3 hybrid scheme is well documented [11,12].
In Figure 4 the relative errors of the computed NO stretching frequencies with respect to the experimental counterparts are plotted for the different families of nitroxide radicals defined above. Even the cheap B3 level of theory performs adequately in the case of acyclic and six-member ring molecules, with comparable values of scaled B3 and B2 harmonic frequencies. However, the agreement between the two levels of theory decreases for nitroxides involving five-member rings (l,m,n,o) with a significant overestimation of ν NO at the B3 level. As a matter of fact, the mean error decreases from 21 (B3) to 6 cm −1 (B2), with the last value representing an excellent achievement when taking into account the variability of the experimental conditions. For all the considered nitroxides the anharmonic corrections are very close (the range being 29-34 cm −1 ) so that a constant correction of 30 cm −1 performs a remarkable job. Alternatively, the computational burden can be reduced by carrying out reduced dimensionality computations in which anharmonic contributions are computed only for normal modes strongly coupled to the NO stretching.
Since a reliable estimate of the coupling is provided by the ratio between the two-mode cubic constants involving the NO mode (K iij , which can be obtained by just two additional Hessian evaluations) and the harmonic frequency of the coupled mode (ω i ), a suitable reduced dimensionality scheme can be selected a priori in terms of a pre-defined threshold. The results collected in Table 5   A final comment is in order about the correlation between the vibrational frequency of the NO stretching and suitable geometrical parameters (e.g., NO bond length or pyramidality of the nitrogen environment). Unfortunately, several attempts to derive sufficiently accurate relationships provided disappointing results possibly because of the delocalization of the NO stretching and the ensuing involvement of several geometrical parameters.

Hyperfine Coupling Constants
The isotropic hyperfine couplings of free radicals are primarily ruled by the shape of the SOMO [9]. For instance, the hcc's of the central atom of π radicals (nitrogen in the case of nitroxides) should vanish for planar structures because of the corresponding lack of any contribution of its s orbitals to the SOMO. The only remaining contribution to a N is spin polarization. As a consequence, as shown by the results collected in Table 6, the a N of planar or nearly planar nitroxides are smaller than those of the pyramidal counterparts. As already mentioned, the correct reproduction of spin polarization is particularly difficult at the DFT level with both high percentages of Hartree-Fock exchange and second order perturbative contributions playing a role to obtain improved results. If this interpretation is correct, the error of DFT results should always decrease with the pyramidality at the radical center due to the lower relative contribution of spin-polarization to the overall a N . The cos 2 (θ) dependence of the error discussed in detail in the section devoted to the template dimethylnitroxide molecule is a direct consequence of this model. On these grounds, two strategies can be employed to improve the computed a N values. In the first one, closely resembling the ONIOM approach [40], the a N computed for the target molecule by a low-level method (here B2) is corrected by the difference of the a N computed for a template molecule (TM, here dimethylnitroxide) at the geometry of the target molecule employing high (here CC) and low level methods. This approach (hereafter referred to as ∆TM) requires a new CC computation of the model system for each target molecule. However, the discussion above suggests an even simpler approach (hereafter referred to as ∆θ) not requiring any additional computation at the CC level. In fact, if the dominant contribution to the error is related to the value of the out-of-plane angle, it is sufficient to fit once for ever the parameters of Equation (1) at the CC level and then use this Equation and the specific value of the out of plane angle for computing the correction to be applied to the B2 value. Table 6. Equilibrium NO bond length (Å), out of plane angle θ (in degrees) and a N (in MHz) computed at the B2 level together with high-level corrections from the template molecule model (∆TM in MHz) and B3 vibrational corrections at 300 K (∆ vib in MHz). In the last column, the correction obtained by Equation (1) (∆θ in MHz) is reported. The column ∆ CNC reports the difference between the CNC angle of the real system and that of the dimethylnitroxide model (a) at the same θ angle. The results collected in Table 6 show that the simplified model provides remarkably accurate results, thus allowing to correct the computed a N of any nitroxide simply in terms of its out of plane angle. In spite of non negligible differences between the CNC angle of the real molecule and that of the dimethylnitroxide model at the same out of plane angle, the discrepancies are always smaller than 0.3 MHz, a value well within the expected error bar of even the most accurate QC models. We have thus at our disposal an effective strategy for computing equilibrium a N values, which can be extended to other density functionals without any additional CC computation, but simply fitting the a and b parameters of Equation (1) for the new functional (by means of cheap computations for the dimethylnitroxide template molecule) to be compared with the already available CC counterparts.
Experimental results actually refer to vibrationally averaged quantities at the temperature of the experiment. These contribution can be effectively obtained in the framework of the GVPT2 model since cubic and semi-diagonal quartic force constants can be computed with sufficient accuracy employing denstity functionals cheaper than double hybrids. In view of previous experience [41], we have employed the B3 model, whose results are also collected in Table 6. It is apparent that the effect of vibrational averaging is particularly significant for some planar or nearly-planar systems, but it cannot be neglected also in other cases in order to obtain quantitative results.
The accuracy of the results collected in Table 6 can be guessed by comparison with the experimental data for the large and particularly demanding INDCO radical p. Summing the contributions of columns 5,6 and 8 we arrive to a N = 25.6 MHz and bulk solvent effects evaluated at the PCM level [22] for the benzene solvent employed in experiments add further 0.5 MHz leading to a final estimate of 26.1 MHz, in remarkable agreement with the experimental value of 25.9 MHz [42].

Concluding Remarks
In the present paper we have proposed and validated a general and robust strategy providing accurate structures, vibrational and magnetic properties of nitroxide radicals by means of hybrid and double-hybrid functionals in conjunction with suitable basis sets. Anharmonicities and vibrational averaging effects on different properties can be taken into account employing the effective generalized vibrational perturbation theory to the second order. The geometrical parameters and vibrational frequencies obtained at this level are remarkably accurate, whereas magnetic properties can be further improved by means of coupled cluster computations performed for a template molecule [8].
However, environmental effects cannot be neglected for studies in condensed phases. For innocent solvents (or more generally when inter-molecular hydrogen bonds are not present) polarizable continuum models can be profitably used and their most effective implementations do not add any significant computational burden [22]. The situation is different for hydrogen-bonding environments since the geometry and magnetic properties of the NO moiety are strongly affected by this kind of interactions. Several studies have been devoted to this problem [4,5,9,10], but a systematic investigation for a large panel of nitroxides and solvents is still lacking. Work is in progress in our laboratory to extend the strategy reported above to those situations by means of effective variationalperturbative approaches resembling those already applied with remarkable success for other spectroscopic properties [43].