## Abstract

We have investigated the structure and nuclear magnetic resonance (NMR) spectroscopic properties of some dihydrogen endofullerene nitroxides by means of density-functional theory (DFT) calculations. Quantum versus classical roto-translational dynamics of H_{2} have been characterized and compared. Geometrical parameters and hyperfine couplings calculated by DFT have been input to the Solomon–Bloembergen equations to predict the enhancement of the NMR longitudinal relaxation of H_{2} due to coupling with the unpaired electron. Estimating the rotational correlation time via computed molecular volumes leads to a fair agreement with experiment for the simplest derivative; the estimate is considerably improved by recourse to the calculation of the diffusion tensor. For the other more flexible congeners, the agreement is less good, which may be due to an insufficient sampling of the conformational space. In all cases, relaxation by Fermi contact and Curie mechanisms is predicted to be negligible.

## 1. Introduction

Guest@host complexes have attracted much attention from diverse scientific perspectives, owing to the unique challenges and opportunities that they offer, ranging from experimental (i.e. synthetic chemistry and spectroscopy) to purely theoretical aspects (geometries and interactions between guest and host). Among these, the dihydrogen endofullerene H_{2}@C_{60} is a veritable ‘gold mine of intriguing aspects for investigation’ [1]: for example, the electronic properties of the C_{60} cage in suitable derivatives may be exploited to trigger ortho–para conversion of the guest H_{2} spin isomers. The paramagnetic relaxation enhancement of molecular hydrogen encapsulated in fullerene C_{60} derivatives (fulleropyrrolidines) bearing a pyrrolidine or piperidine oxyl functionality (endofullerene nitroxides) has been investigated [2]. A series of such nitroxides that have been investigated is displayed in figure 1. In these compounds, the paramagnetic contribution to the proton nuclear magnetic resonance (NMR) longitudinal relaxation time (*T*_{1,p}) of encapsulated H_{2} was measured and found to approximately follow a 1/*r*^{6} law, *r* being the distance between the nitroxide functionality and the centre of the C_{60} cage, where the H_{2} molecule was supposed to be constrained. This result suggests that the dipole–dipole mechanism is responsible for relaxation, provided that the nucleus–electron spin–spin distance (*r*_{SI}) can be approximated to *r*, as depicted in figure 1*b*.

Despite the overall consistency of the results, some issues have remained open: (i) the magnitude of the values of *r* were estimated only by molecular mechanics calculations, and hence are not very accurate; (ii) paramagnetic relaxation was assumed to take place only through nuclear–electronic spin dipolar coupling (i.e. Fermi contact and Curie relaxation were not considered); and (iii) correlation times for the rotational motion of the whole molecule were estimated with a very simplified model. In this paper, we wish to address some of the above issues for the molecules of figure 1*a*.

First-principles electronic structure calculations, notably those based on density-functional theory (DFT), are capable of predicting the parameters related to electron paramagnetic resonance (EPR) and NMR spectroscopies [3,4,5,6]. Whereas the calculation of the NMR properties of diamagnetic molecules is widely applied [3–7], the corresponding computations for paramagnetic species are comparatively less developed. Nevertheless, it has been shown that DFT methods can be used to predict NMR and EPR parameters for a variety of paramagnetics [8–20].

In fact, a bigger stumbling block specific to modelling H_{2} (or other small molecules) encapsulated in C_{60} or its derivatives is related to the quantum nature of its internal dynamics. Bačić and co-workers [21] have extensively investigated the problem up to five-dimensional roto-translational (RT) quantum calculations, adopting an empirical interaction potential for dispersion. Even before H_{2}@C_{60} was synthesized, Cross had studied the rotation of encapsulated H_{2} with a potential based on a Hartree–Fock wave function empirically corrected for dispersion [22]. Those efforts (see [23] for a review) enabled a full description of the quantum RT energy levels, thus describing the motion of H_{2} inside the cage—also in connection with inelastic neutron scattering (INS) experiments [24]. However, apart from the simplified potential adopted, these approaches do not provide the electronic structure of the H_{2}–C_{60} ensemble, which is obviously essential for predicting and modelling NMR/EPR parameters and relaxation. Such calculations are generally carried out with high-quality quantum-chemical methods on fixed geometrical arrangements, which are obtained by gradient optimization or time snapshots of molecular-dynamics runs. In the latter case, interaction potentials may be classical (force-field based) or quantum (Car–Parrinello and variants thereof), but RT motion is always treated classically. As mentioned earlier, this approach is questionable for H_{2}@C_{60}. A case study concerning the related Xe@C_{60} system was provided by Straka *et al.* [25], where the NMR of encapsulated ^{129}Xe was investigated.

We are then left with the problem of how to deal with the electronic structure of a system whose internal dynamics is dominated by quantum effects. It has even been argued that any attempt at treating H_{2}@C_{60} as a rigid entity with fixed atom positions is ‘inherently wrong’ [26]. On the other hand, it is conceivable that the population of RT energy levels selects a region of space where H_{2} can be *operationally* located, i.e. ascertain whether certain H_{2}–C_{60} fixed arrangements can be used as ‘placeholders’ to investigate the electronic structure of the dynamic ensemble; this is the avenue we pursue here.

The goals of this study are manifold, namely: (i) elucidate the energetics of arrangement of the encapsulated hydrogen molecule with respect to the fulleropyrrolidine structure, especially in relation to H_{2} dynamics; (ii) understand the effect of such structural features on the NMR-relevant properties; (iii) improve the estimated distance, *r*, between the radical and C_{60} centres, and justify the previous assumption that the hydrogen nuclei can be localized at the cage centre for the purpose of modelling the paramagnetic relaxation enhancement; and (iv) assess the importance and contributions, if any, of the Fermi contact and Curie relaxation. In this context, an accurate estimate of the rotational correlation times of the whole molecule will turn out to be of critical importance.

The results are divided into two main sections: §2*a* concerns the dynamics of H_{2} in C_{60} as a model for that of H_{2} in fulleropyrrolidines, and §2*b* concerns the electronic structure and NMR/EPR properties of H_{2} within compounds **1–5** of figure 1*a*.

## 2. Results and discussion

### (a) Model for longitudinal nuclear spin relaxation

#### (i) Sketch of the formal approach

The paramagnetic contribution to the longitudinal relaxation rate (*R*_{1,p}=1/*T*_{1,p} as defined in [2]) stems, to different extents, from the modulation of dipolar and Curie interactions, and possibly also from the Fermi contact interaction.

The formal approach yields the NMR signal evolution through the solution of the phenomenological stochastic Liouville equation (SLE) [27], including the effect of the radiofrequency pulses in the inversion–recovery pulse scheme, for the evolution of the spin density matrix under modulation of the spin–spin interactions. Along this line, one should start from the representation of the density matrix operator of the whole spin system (nuclear–electronic) over the quantum states, which are, *a priori*, spin states and intra-cage RT states of H_{2}. Spin relaxation comes from stochastic modulation (of the spin Hamiltonian), which enters at two levels. First, fluctuations of the classical degrees of freedom are accounted for in the SLE through a multidimensional Fokker–Planck operator [28] for the intramolecular conformational dynamics and the tumbling of the whole molecule in the isotropic liquid phase (diffusive processes); the modelling of such dynamics requires the detailed knowledge of the internal energetics, conformation-dependent hydrodynamics and consideration of recoil effects [29]. This would cover the whole range of situations from the rigid limit of distinct conformers, up to the opposite limit of very fast structural fluctuations within potential wells plus conformational transitions; while compound **1** falls in the former category due to its rigidity, the latter situation is likely to apply for all molecules **2**–**5**, albeit to different extents. A further source of stochastic modulation comes from transitions among the intra-cage RT levels, which can be accounted for in the SLE by means of a relaxation matrix inserted *ad hoc*. We stress that algebraic elaboration and numerical solution of the full form of the SLE is not only quite demanding, but its usefulness is questionable as well, since many of the input ingredients are unknown. In particular, hereafter we shall briefly discuss the problem concerning the RT transition matrix.

#### (ii) Intra-cage roto-translational quantum states and dynamics

Although the finite extension of the linear H–H molecule breaks the icosahedral symmetry of the C_{60} cage, RT coupling can be neglected to a first degree of approximation. Recent INS experiments [24] have shown that H–H in C_{60} can be approximated to a free quantum rotator. The energies of the quantum states are expressed by , *J*≥0 and −*J*≤*M*_{J}≤+*J* being the principal and projection angular momentum quantum numbers, respectively. The rotational constant is *B*≈0.285 *k*_{B}*T* at 298 K (from the energy separation of 14.6 meV between the *J*=0 and *J*=1 levels) [24]. For completeness, both spin allotropes of molecular hydrogen are considered here, namely *o*-H_{2} (ortho) and *p*-H_{2} (para), although only *o*-H_{2} is detectable by NMR. We recall that for *o*-H_{2} and *p*-H_{2} only levels with odd or even *J* values are allowed, respectively: with this rule, the essential set of rotational quantum states can be assessed by looking at the rate of convergence of the Boltzmann summation in approaching the value of the quantum rotational partition function . This is illustrated in figure 2: three states suffice for the ortho form (the three degenerate states with *J*=1) and six states are required for the para form (the single state with *J*=0 and the five degenerate states with *J*=2).

To tackle the three-dimensional translational problem for the H_{2} centre of mass, we have first generated one-dimensional profiles of the potential energy versus the displacement from the cage centre for two limiting orientations of the H_{2} molecule: collinear and orthogonal to the displacement direction. Potential energies were calculated with the PBE0 density functional, including dispersion (PBE0-D3/TZ2P level; see §§2*b*(i) and 4 for a full description).

The outcome is shown in figure 3*b*, where continuous lines represent a best fit with polynomials of order up to six. The similarity between the profiles supports the view of little RT coupling, at least at the level of description adopted here. The simple average of these profiles is taken as a rough choice for the effective radial potential of spherical symmetry to be used within the approximation of RT decoupling.

The three-dimensional stationary Schrödinger equation has been solved in spherical coordinates for such an average potential. The usual angular–radial decoupling between azimuthal/polar angles (*θ*_{c},*φ*_{c}) and distance *r*_{c} of the centre of mass of H_{2} from the cage centre yields eigenfunctions factorized as *ψ*_{l,m,n}(*θ*_{c},*φ*_{c},*r*_{c})=*Y* _{l,m}(*θ*_{c},*φ*_{c})*R*_{n,l}(*r*_{c}), where *Y* _{l,m}(*θ*_{c},*φ*_{c}) are spherical harmonics with *l*, *m* the principal and projection quantum numbers, and *R*_{n,l}(*r*_{c}) is the radial component, *n* being a progression index. Following the standard route of treatment of central-field quantum problems, solution of the set of *l*-dependent radial equations has been achieved by means of a finite-difference scheme with boundary placed at 1.75 Å(internal radius of the C_{60} cage estimated from the van der Waals radius of carbon atoms); ranks up to *l*_{max}=10 and 100 radial intervals yielded full convergence on any property calculated here. The resulting quantum states are (2*l*+1)-fold degenerate on *m*, with energies , which have been ordered by increasing magnitude. The first terms of the series are presented in figure 3*c*. As for the rotational problem, the convergence of the Boltzmann summation to the asymptotic value (the translational quantum partition function) is taken as an indicator to select the minimal set of translational states to be handled in the SLE. The outcome is shown in figure 3*d*, and suggests that an extended pool of states, up to several tens, should be taken into account. The much smoother growth of the profile in figure 3*d* with respect to those for the rotational case is not surprising, since for the translational case the states appear in groups with almost constant spacing along their progression, and the degree of (quasi-)degeneracy grows more rapidly as the energy increases. The overall behaviour is qualitatively similar to the progression of RT states for H_{2}@C_{60} calculated from five-dimensional quantum dynamics [21]. Many of these highly degenerate states are individually scarcely populated by thermal fluctuations due to their high energy, but their collection plays an important role in the overall statistical analysis.

This analysis indicates that, at room temperature (298 K), the description of the RT intra-cage dynamics of H_{2} involves a large number of quantum states (hundreds of RT states when taking the direct products of translational and rotational states). Moreover, appreciable energy spacing of the order of *k*_{B}*T* places the encapsulated H_{2} right at the border between the quantum and classical descriptions. Nevertheless, for practical purposes, the proper quantum equilibrium distribution at room temperature can be replaced by the classical expression, which is much easier to handle in the evaluation of the quantities required for the present work. Thus, as detailed in appendix A, the quantum and classical translational distributions are barely distinguishable at 298 K (our reference operating temperature), while the quantum and classical rotational distributions are identical. This finding has some important implications: despite the fundamentally quantum nature of encapsulated H_{2} dynamics, at least at room temperature the properties of such systems can be understood (and reliably dealt with) by means of a classical description. It should be remarked, however, that such an approximation breaks down at lower temperatures (the case of 50 K is explicitly shown in appendix A).

#### (iii) Simplified approach to relaxation

Specification of the relaxation matrix to be input in the SLE requires detailed modelling (currently not available) of the interaction of such an open quantum system with the many classical or near-classical environmental modes (i.e. cage vibrations, structure and local solvation of the fulleropyrrolidine) acting as thermal bath. Given the large number of involved RT states, an approximate classical description of the RT motion can be worked out in terms of an ‘effective’ semi-inertial dynamics of H_{2} in a low-friction environment. RT motion can then be treated by means of a probabilistic Kramers–Klein equation [28,30] in which, however, the input friction coefficients within the C_{60} cavity have to be determined *a priori* (only a detailed analysis of wave-packet evolution for the RT quantum dynamics can provide such information). This problem, together with the difficulty of fully accounting for the conformational dynamics, calls for a simplified approach. Therefore, regardless of the details of RT and conformational dynamics, we shall apply a simplification based on the assumption of time-scale separation: both conformational motions (if present) and RT dynamics of H_{2} are much faster than the tumbling of the molecule as a whole; moreover, hydrogen dynamics is reasonably much faster than any conformational process. Accordingly, the spatial terms in the spin Hamiltonian are progressively modulated in such a way that effective ‘residual’ Hamiltonians are left for modulation by slower dynamics. Such a dynamical decoupling allows for the paramagnetic contribution to the longitudinal relaxation rate to be expressed as the sum of independent terms
2.1where ‘rot’ stands for the tumbling process of the whole system. The single terms are evaluated separately by means of suitable models for the corresponding motions. For example, can be expressed in the frame of the well-known model-free (unspecific) approach of Lipari & Szabo [31,32] if the internal dynamics is assimilated to a ‘wobbling in a cone’ of a small mobile fragment with a single dominant correlation time. For , as indicated in §2*a*(ii), a fast semi-inertial RT motion of H_{2} could be a likely choice, although a justification has to be provided on sound grounds. The term due to the molecular tumbling, on the contrary, is well characterized in the framework of Solomon–Bloembergen equations as shown in §2*a*(iv).

A step further is to assume that, to a first approximation, the last two terms in equation (2.1) are negligible compared with . This holds true, generally speaking, as long as the conformational and intra-cage motions are fast and ‘localized’, in the sense of being poorly efficient in modulating the spatial part of the spin Hamiltonians.^{1} Note that the assumption of negligible can be checked *a posteriori* on compound **1**, for which conformational dynamics is absent. We anticipate that such a contribution is indeed negligible for compound **1**, as shown in the following; the same must hold true for compounds **2–5**. We shall now focus on the supposed leading term in equation (2.1) due to molecular tumbling as the main process in driving nuclear relaxation.

#### (iv) Tumbling contribution

Theory for the relaxation rate originated by isotropic molecular tumbling is well established. It is expressed as a sum of individual contributions due to dipolar, Fermi contact (FC) and Curie interactions,
2.2where each term is given by the following Solomon–Bloembergen equations [33]:
2.3
2.4
and
2.5Here *γ*_{I} is the magnetogyric ratio of the resonant nucleus, *g*_{e} the isotropic *g*-factor, *S* the spin quantum number, *A* the isotropic hyperfine coupling constant (in hertz), *ω*_{I} and *ω*_{S} the nuclear and electronic Larmor frequencies (with |*ω*_{I}|≪|*ω*_{S}|), *τ*_{S} and *τ*_{r} the electronic relaxation time and the rotational correlation time, respectively, and *τ*_{c} the overall correlation time, which, neglecting any exchange mechanism, can be expressed as
2.6Since *τ*_{S} for fulleropyrrolidine nitroxides is *ca* 10^{−6} s [33–35] and *τ*_{r} lies between 10^{−10} and 10^{−9} s (see below), then *τ*_{c}≈*τ*_{r}. For molecules **1**–**5**, *S*=1/2 and, at a given magnetic field, *ω*_{I} and *ω*_{S} are also known; herein we will adopt a proton Larmor frequency *ω*_{I}/2*π*=500 MHz [2]. Finally, we recall that *r*_{SI} in equations (2.3) and (2.5) accounts for the distance between the electron and the resonant nucleus.

According to the time-scale separation invoked above, molecular tumbling modulates the ‘residual’ Hamiltonians left by the action of faster motions; hence the average 〈*r*^{−6}_{SI}〉 is meant to be taken over the intra-cage RT equilibrium distribution of the H_{2}, and over the conformational equilibrium distribution. However, since Boltzmann averaging over the conformational space is highly demanding on computational grounds, only the most stable conformer will be taken into account. On the contrary, the RT average is made explicitly, adopting a classical, instead of quantum, intra-cage RT equilibrium distribution, as mentioned earlier (see appendix A). The angular distribution of the H–H axis is taken to be isotropic, while the Boltzmann distribution of the centre of mass is constructed by using the effective potential discussed earlier. Recalling that *r* denotes the distance of the electron spin from the centre of C_{60}, the average can be expressed as
2.7where *f*_{T} is a dimensionless parameter dependent on *r*. Calculations of *f*_{T} versus *r* have been done by numerical averaging; a value 0.8 Åhas been assumed as the distance between the protons in the hydrogen molecule. It turns out that *f*_{T} decays monotonically as *r* increases, and it differs very slightly from 1 (*f*_{T}≈1.10 at *r*=4 Å and 1.01 at *r*=15 Å). For compounds **1**–**5**, where *r* is larger than 6 Å, one can safely replace 〈*r*^{−6}_{SI}〉 with *r*^{−6} with less than 4% error.

The essential conclusion is that averaging over the intra-cage distribution allows one to proceed *as if each resonant nucleus was placed at the centre of C*_{60}. This outcome justifies the previously adopted choice of taking *r* as the electron spin–nuclear spin distance in the Solomon–Bloembergen equations [2].

In conclusion, to evaluate , the unknown parameters for each molecule of the series are *g*_{e}, *A* and *r*, along with the rotational correlation times *τ*_{c} that depend to a large extent on the bulk hydrodynamic properties of the molecule in the solution. We stress again that, instead of performing a highly demanding Boltzmann averaging over the conformational space, the effective distance *r* and the rotational correlation time *τ*_{c} are evaluated for the most stable conformer.

### (b) Electronic structure calculations

#### (i) Computational issues

Modelling the electronic structure of endofullerene nitroxides (open-shell systems with more than 60 non-hydrogen atoms) raises several technical issues, because (i) the computational size is substantial even for the smallest congener (**1**) and (ii) the energetics of interaction of encapsulated H_{2} with the fullerene cage is dominated by weak dispersion forces, which are not correctly modelled by standard density functionals. The same issues hold also for the model H_{2}@C_{60} system examined before (see also [36,37]). The need to fulfil both requirements led us to choose a calculation protocol based on Slater basis sets with a GGA (BP86) or hybrid (PBE0) functional, corrected for dispersion energies (see §4).

A further issue to deal with concerns the orientation of the H_{2} molecule with respect to the fulleropyrrolidine ring bearing the unpaired electron. Since the potential energy surface (PES) for motions of H_{2} close to the centre of C_{60} is quite flat (figure 3), it is likely that the forces acting on H_{2} are too small to allow for a reliable convergence of geometry optimization also for fulleropyrrolidines. Therefore, in order to address this sensitive issue (a change in orientation entails a change in *r*_{SI} for each nucleus), we investigated endofullerenes **1**–**5** with six starting geometries where H_{2} is oriented along the *x*-, *y*-, *z*-axes (the *z*-axis pointing towards the NO group), as well as along the axes defined by the bisector of the *xy*, *xz* and *yz* planes (*xy*, *xz* and *yz* ‘axes’ for short; figure 1). In all cases, the centre of mass of H_{2} was placed at the centre of the fullerene cavity.

#### (ii) Compound **1**

As a first step, the BP86 functional was selected, because it leads to much faster calculations than with hybrids, and because the *epr* module in ADF (the Amsterdam Density Functional program) allows for the calculation of *g*-factors only with GGA functionals. All optimized structures retained their initial arrangement of H_{2} relative to the pyrrolidine ring, and the PES was indeed extremely flat, as indicated by the difficulty in achieving geometry convergence even with tight criteria (10^{−4} atomic units (Hartree) in energy). All arrangements featured quite similar energies within 0.1 kcal mol^{−1} of each other and can be considered isoenergetic for all purposes (although the *x* orientation was consistently more favoured). All *g*-factors were also essentially identical, with an average value of *g*_{e}=2.0067, typical of nitroxides [38].^{2}

Since hybrid functionals are generally regarded as superior to GGA ones, we benchmarked the energies obtained at the PBE0-D3/TZ2P level with both the BP86-D3 and PBE0-D3 geometries (the latter series of calculations being quite demanding in the optimization step). The results are qualitatively still similar, with extremely small energy differences with orientation. The PBE0-D3//BP86-D3 calculations exhibit the largest range, which however remains below 0.5 kcal mol^{−1}, i.e. well within the expected error of these methods, also considering the difficulties in geometry convergence. The relative energies of all six orientations are collected in figure 4.

Hyperfine couplings were averaged over the six H_{2} orientations according to Boltzmann populations. This is only a rather rough way of accounting for the actual orientational distribution of the H_{2}; however, the results did not differ much from an arithmetic average because calculated populations were obviously quite similar.

As expected, the relative magnitude of *A* for each ^{1}H spin differs according to the orientation with respect to the nitroxide moiety: thus, in the *x* arrangement *A* values are equal, whereas along *y*, *z* and their combinations the values of *A* become different in magnitude and sign (figure 5; and see footnote 2). However, the absolute values do not exceed *ca* 5 kHz in the *yz* orientation; averaging according to the populations and both hydrogen atoms, we obtain *A* values between 0.2 and 6 kHz according to the method used. Hereafter, the average value of *A*=1.6 kHz will be adopted for the estimation and comparison of relaxation rates of **1**–**5**; indeed, any value within that range would lead to the prediction that the hyperfine mechanism is negligible for **1–5** (see below).

A similar conclusion was reached by comparing *T*_{1} with an estimate of *T*_{2} for compound **1** [40]. Therein, it was also concluded that a scalar contribution to *T*_{2} would have been detectable for a value of *A* greater than *ca* 3 kHz. This is consistent with the values of *A* found in this study.

The frontier molecular orbitals of **1** are displayed in figure 6 (in the *x* orientation by way of example) and are spatially very different. The highest, singly occupied molecular orbital (HOMO) is localized on the nitroxide, as expected, whereas the lowest unoccupied molecular orbital (LUMO) essentially spans only the fullerene cage. The HOMO–LUMO gap (*ΔE*) is heavily dependent on the level of theory: with BP86 *ΔE*=0.0161 Hartree, whereas with PBE0 it is much larger (*ca* 0.09 Hartree, slightly depending on the geometry used). This feature caused severe self-consistent field (SCF) convergence problems with **2**–**5**, as reported below. The spin density isosurface of **1** closely resembles that of typical nitroxides, being localized essentially on the N and O atoms [12].

#### (iii) Compounds **2**–**5**

Apart from the larger molecular size, the most prominent feature of the remaining endofullerenes is the steadily decreasing HOMO–LUMO gap computed with the BP86 functional: *ΔE* varies from 0.0161 Hartree (**1**) to a vanishing 4×10^{−5} Hartree (**5**). Such small gaps, joined with the spatial structure of the HOMO and LUMO seen above, led to severe SCF convergence problems. Algorithms designed for such problems allowed in some cases for successful convergence, but with huge computational expense. Most importantly, the very small *ΔE* was found to be an artefact of BP86 (indeed, of most GGA functionals; data not reported). Analogous calculations with PBE0 (or B3LYP) led to *ΔE* values of 0.06–0.10 Hartree, i.e. similar to those for **1** and rather constant along the series. While the obvious solution would be to run all calculations with PBE0, geometry optimization would be too expensive, mainly owing to the very large number of optimization steps required. As a workaround, we optimized the geometries with BP86 and level shifting (see §4), and computed spectroscopic parameters with PBE0 at that geometry. Optimized geometries featured the conformations sketched in figure 1, with the cyclohexane rings of the tetramethylpiperidine moiety in the chair conformation that places the methyl groups farthest from the fullerene surface (however, see below). Derivative **5** may also exhibit several conformations deriving from rotation about the C−C bond connecting the pyrrolidine and the phenyl ring. Some other such conformations were tested, and the one displayed in figure 1 was found to be the most stable.

As mentioned, we could not compute *g*_{e} with this method; however, in the calculations where BP86 converged, the values were (not surprisingly, in view of the close similarity) essentially identical to those of **1**; therefore, we adopted that value (2.0067) through the series.

#### (iv) Estimating the paramagnetic relaxation rate of compound **1**

The distance *r* was set to 6.72 Å(from the BP86 geometry for consistency with further calculations; see below). The molecular volume *V* was estimated at 4623 Å^{3} and the viscosity was set at *η*=1.322 cP (the value at room temperature for *o*-dichlorobenzene, the solvent employed for the measurements [2]). Also, we assumed *τ*_{S}=10^{−6} s, as typical for nitroxides. Applying the Debye–Stokes–Einstein model (as in [2]), one obtains *τ*_{r}=1.5 ns for the overall rotation of the host cage. Finally, setting *g*_{e} to the average value of 2.007, we obtain , and . Quite obviously, relaxation takes place exclusively by the dipolar pathway, so that *T*_{1,p}=19 ms (to be compared with the experimental value of 13.5 ms).

The estimate of *τ*_{r} was improved by calculating the full diffusion tensor (see §4), whereby a refined value of *τ*_{r}=0.69 ns (with stick boundary conditions at the fulleropyrrolidine surface; see below) was obtained. Using the same parameters as above, we now have , that is *T*_{1,p}=10.2 ms, in much better agreement with experiment. Hereafter, correlation times calculated from the diffusion tensor will be used.

#### (v) Paramagnetic relaxation rates of **1**–**5**

In the calculation of *τ*_{r}, two limiting boundary conditions for the shear regimes of the surrounding fluid at the molecular surface can be adopted, namely ‘stick’ and ‘slip’. The standard choice is stick conditions, although slip has been suggested to be appropriate for almost spherical molecules such as fullerenes [41]. A trial slip calculation for **1**, the closest to a spherical object, leads however to a shorter *T*_{1,p} (8.3 rather than 10.2 ms), thus making the agreement with experiment worse. Since there is no cogent rationale for adopting slip conditions, only results from stick conditions will be presented. In table 1, we report estimates of *τ*_{r} and *T*_{1,p} for stick conditions only; in all instances, relaxation turns out to be essentially dipolar.

With reference to the experimental values, the calculated relaxation times are underestimated by *ca* 25% for **1**, overestimated by 50% for **2** and **3**, and even more for **4** and **5**. The corresponding overestimation of the relaxation rate for **1** is anomalous, since even an incomplete but well-parametrized model should lead to an underestimation (contributions are additive). However, the overall agreement for **1** is good, considering the approximations involved in the calculation and the errors inherent in the experimental determination ( is obtained by subtraction of the relaxation rate of the corresponding diamagnetic hydroxylamine, i.e. *R*_{1,p}=*R*_{1}−*R*_{1,dia}) [2].

The dependence of the relaxation rate on *τ*_{c} and *r*_{SI} (equation (2.3)) is shown in figure 7. The relevant *τ*_{c} values lie in the region of maximum steepness for equation (2.3), where the relative error *ΔR*_{1}/*R*_{1} has the highest sensitivity to *τ*_{c}.

The good agreement found for the most rigid compound **1** supports the conclusion that correlation times are estimated correctly using stick boundary conditions. Moreover, the outcome supports the above simplification of neglecting the dynamical effects of the H_{2} rotational motion on *R*_{1,p}. In fact, the good accord (and even the slight underestimation of the rate) points to the statement that the excluded term in equation (2.1) is indeed immaterial. This *a posteriori* check validates our methodological assumption, which then applies to all other compounds of the series.

Conversely, for **2**–**5** the agreement is worse. A rationale for the degrading performance going down in the series (up to a factor of 2.3 in terms of *T*_{1,p} for **5**) goes as follows. The discrepancies cannot be attributed to a poor estimate of *τ*_{r} (and hence of *τ*_{c}) under stick boundary conditions, in view of the good agreement for compound **1** seen above. The discrepancy for **2**–**5** may then be attributed to their higher conformational flexibility, which has not been taken into account so far. First, an improvement of the estimate for would be obtained with *r* resulting from a proper average over the conformer's distribution; second, the other excluded term in equation (2.1), that is, , may play some role and needs to be inspected.

In order to discuss conformational effects, we begin by focusing on compound **2**, which is still relatively rigid but has another conformer **2**′ where the piperidine ring has undergone chair inversion (figure 8). This arrangement features a closer approach of the methyl groups to the fullerene surface and is accordingly less stable by *ca* 5 kcal mol^{−1} than **2**. This implies that such a configuration has a negligible weight in the thermal average (of the order of 10^{−4} compared with **2**) and hence it is essentially not sampled by thermal fluctuations. On the other hand, already for this rigid compound we have a mismatch of *ca* 50% for the relaxation rate (time), which cannot be reasonably attributed to small-amplitude and very fast librations of the nitroxide moiety.

To get a deeper insight into this feature, we can consider the sources of error in the estimate of *R*_{1,p}. From equation (2.3), the relative error in *R*_{1,p} is given by *ΔR*_{1,p}/*R*_{1,p}≈6*Δr*/*r*+*Δτ*_{c}/*τ*_{c}, whereby it is apparent that the relative error in the distance *r* contributes more than the relative error in the correlation time *τ*_{c}. If one assumes that *Δτ*_{c}/*τ*_{c} is negligible, the relative error *ΔR*_{1,p}/*R*_{1,p}≈50% found for **2** can be accounted for by an 8% error in the estimate of *r*. Given that the calculated distance *r* is about 9 Å, a misplacement of *ca* 0.7 Åcan account for the observed discrepancy. While this uncertainty may seem large, one should keep in mind that *r* is not defined rigorously. In keeping with our previous work [12], we have taken *r* as the distance between the centres of the internuclear vectors in the H–H molecule (i.e. the centre of the C_{60} cage) and in the N–O group, which proved to be a reasonable working solution. However, this choice does not account for the delocalization of the electron spin. This turns out to be critical in the present case, where relatively small discrepancies are found; alternative approaches, not based on the point-dipole approximation [10], may prove helpful. Therefore, we can conclude that the predicted relaxation rate for **2** is consistent with experiment within the understood sources of error.

Compounds **3**–**5** exhibit even larger discrepancies. For these flexible molecules, in addition to the sources of error discussed for **2**, the contribution from conformational dynamics would have to be accounted for. A detailed evaluation of such internal dynamics is, however, beyond the scope of the present work.

## 3. Conclusions

The calculation of the molecular and spectroscopic properties of H_{2} endofullerenes poses several fundamental problems in addition to known computational issues. These problems are related to the quantum nature of the translational and rotational dynamics of the encapsulated molecule, which in principle clash against the need for fixed rigid arrangements required for accurate prediction of the electronic structure. Nevertheless, we have shown that (at least near room temperature) the conventional approach of dealing with a finite set of fixed hydrogen/cage arrangements is an effective working model. Thus, we have computed, by means of state-of-the-art DFT methods and hydrodynamic models, the parameters (geometrical, electronic and dynamic) required for the estimation of the paramagnetic relaxation rate of molecular hydrogen encapsulated in fullerene nitroxide derivatives. Relaxation is by far dominated by the dipolar pathway. The calculated estimates of the distance between hydrogen and the nitroxide functionality (*r*) and the rotational correlation time (*τ*_{r}) have allowed for a very good agreement with experiment as regards the relaxation rate of the conformationally rigid derivative **1**. For the remaining, more flexible, derivatives, the agreement is worse; a detailed evaluation of the internal dynamics is necessary and will be the subject of future studies. However, apart from the specific instance presented herein, this study opens the way for extensive computational investigations of the intriguing chemistry and physics of endofullerenes.

## 4. Computational details

### (a) Electronic structure calculations

All calculations were carried out with the ADF 2010 suite of programmes [42] with the BP86 [43,44] or PBE0 [45] functionals. Dispersion energies were modelled via Grimme D3 correction [46]. A triple-zeta, twice-polarized basis set made up of Slater functions (TZ2P) was employed. The computational levels (spin-unrestricted) are then denoted as BP86-D3/TZ2P (or PBE0-D3). Hyperfine couplings (*A*-tensors) were computed at the optimized geometry; *g*_{e} factors were computed with the ADF *epr* module [42,47]; however, this has been implemented only with generalized gradient approximation (GGA) functionals such as BP86. For reasons detailed in the text, geometry optimization was in most cases carried out only at the BP86-D3/TZ2P level, applying a level shifting of 0.03 Hartree to ensure SCF convergence. Calculation of energies and *A*-tensors was then performed at the PBE0-D3/TZ2P level.

Each calculation provided the geometry, energy, *A*-tensors of encapsulated H_{2} and *g*_{e}, as well as the molecular volume. The distance parameter *r* was determined as the distance between the centre of the C_{60} cage and the midpoint of the N−O bond for the optimized structure [12,13].

The Cartesian coordinates, energies and spectroscopic properties of all the species investigated are available from the authors upon request.

### (b) Calculation of the correlation time

Rotational correlation times were first estimated through the rough relationship *τ*_{r}=*ηV*/*k*_{B}*T*, computing the molecular volume *V* from the solvent-excluded surface. More accurate evaluation of *τ*_{r} was done through the calculation of the full RT diffusion tensor for the fullerene cage at atomistic level. The molecules were built from overlapping atomic beads, and the RT friction tensor was evaluated with inclusion of hydrodynamic interactions. Stick conditions for relative flow at the molecule–solvent boundary was assumed, except where noted in the text. The principal axes and components (, and ) of the rotational block of the diffusion tensor were then evaluated according to Carrasco & de la Torre [48]. The rotational diffusion tensor was found to be approximately axially symmetric, i.e. and , with the longitudinal principal axis approximately collinear with the vector connecting the centre of the C_{60} cage with the nitroxide group (*z*-axis in figure 1). This allowed us to express the correlation time as .

## Funding statement

Financial support from the University of Padova (PRAT CPDA103095) is gratefully acknowledged. The authors at Columbia thank the NSF for its generous support of this research through grant CHE 1111392. Calculations were run on the Linux cluster of the Laboratorio Interdipartimentale di Chimica Computazionale (LICC) of the University of Padova.

## Acknowledgements

Conversations with Prof. Z. Bačić are gratefully acknowledged.

## Appendix A. Quantum and classical RT equilibrium distributions of H_{2} inside the C_{60} cage

Under the assumption of RT decoupling in the spherical potential, the Boltzmann equilibrium distribution of H_{2} for positional and orientational coordinates (of its centre of mass with respect to the centre of C_{60}, **r**_{c}, and spherical angles *θ*, *φ* of the H–H axis) is factorized as follows:
A1where the two components are normalized separately.

Here, we demonstrate that, for practical purposes, the proper quantum equilibrium distribution can be replaced, at room temperature (298 K), by the classical expression, which is much easier to handle in order to evaluate the averages required in the present work (namely, the averages 〈*r*^{−6}_{SI}〉 that enter equations (2.3) and (2.5)).

First, we consider the translational term *ρ*_{T}(**r**_{c}). Owing to the angular–radial decoupling, we can focus separately on the radial and angular components, *ρ*_{T,rad}(*r*_{c}) and *ρ*_{T,ang}(*θ*_{c},*φ*_{c}), respectively. Let us consider, for the moment, only the radial part normalized as (with *R*_{max} giving the physical boundary, that is, the internal radius of the cavity). The proper quantum distribution is expressed by a weighted sum of squared modulus of the radial components *R*_{n,l}(*r*_{c}) of the eigenfunctions resulting from the numerical solution of the Schrödinger equation (see §2*a*(ii) for details and notation). Namely,
A2with the temperature-dependent weights *w*_{n,l} given by the Boltzmann factors
A3where the degeneracy on the projection quantum number *m* has been explicitly accounted for. The classical counterpart is evaluated, using the average radial potential (see §2*a*(ii)), as the Boltzmann distribution
A4Figure 9 reveals that the calculated profiles from equations (A2) and (A4) are almost identical at 298 K, which is the operating temperature for ref. [2]. The inset shows that only at very low temperatures (here we show the case 50 K) are the quantum and the classical profiles clearly distinguishable.

Let us now consider the rotational component *ρ*_{R}(*θ*,*φ*) in equation (A1), and demonstrate the equivalence of the quantum and classical expressions. The same arguments then identically apply to the analogous problem, still pending, concerning the angular component *ρ*_{T,ang}(*θ*_{c},*φ*_{c}) of the translational term *ρ*_{T}(**r**_{c}). *A priori* one can state that, under the assumption that H–H behaves as a free rotator (see §2*a*), such a distribution is constant, that is
A5at *any* temperature and regardless of whether the quantum expression or the classical approximation is adopted. Intuitively, this comes from the fact that any orientation of the H–H axis *must be a priori* equally probable. On more formal grounds, the quantum expression (the analogue of equation (A2) for the rotational problem)
A6always gives 1/(4*π*) for any values of angles *θ*, *φ*. In fact, *ψ*_{J,MJ}(*θ*,*φ*) in equation (A6) are the eigenfunctions of the quantum free rotator, that is, the normalized spherical harmonics of rank *J*≥0 (the restriction to even values has to be inserted to specifically treat the *o*-H_{2} form important in the NMR context) and −*J*≤*M*_{J}≤+*J*. The weight factors are evaluated as in equation (A3): note that *w*_{J,MJ} are independent of *M*_{J}, since for the free rotator the rotational states −*J*≤*M*_{J}≤+*J* are degenerate. Hence equation (A6) can be rewritten as
A7where the new factors are with by normalization. Since the internal sum over *M*_{J} yields from the spherical harmonics addition theorem [49], equation (A5) follows immediately. From the classical point of view, equation (A5) holds trivially for the unbiased rotator.

Summarizing, at room temperature, we introduce no appreciable errors in replacing the proper quantum expressions with the classical approximation for the RT equilibrium distributions of H_{2}.

## Footnotes

One contribution of 13 to a Theo Murphy Meeting Issue ‘Nanolaboratories: physics and chemistry of small-molecule endofullerenes’.

↵1 For example one can state that results in a weighted sum of terms like those within square brackets of equation (2.3), each term being related to an independent fluctuation mode. For each mode, the correlation time is the inverse of the characteristic relaxation rate associated with it (strictly speaking, an eigenvalue of the diffusion operator for the conformational dynamics, which can be related to the kinetic constants for transitions between stable conformers), and the associated weight factor is proportional to the squared projections of the modulated dipolar component of the spin Hamiltonian over such a mode. Extremely fast conformational transitions (very short correlation times) and/or small-amplitude displacements (small weight factors) point towards being negligible with respect to .

↵2 The difference in

*A*between the two protons for some orientations may also be described as a field gradient across the H_{2}molecule of the sort required for conversion between the ortho and para forms by the Wigner mechanism. The calculations for compound**1**suggest that a hyperfine gradient may make some contribution, in addition to the dipolar interaction, to the conversion rate of**1**, which has so far proved too fast to measure by the methods currently available [39].

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.