Fullerene C 70 characterization by 13 C NMR and the importance of the solvent and dynamics in spectral simulations †

The nuclear magnetic resonance (NMR) spectroscopy combined with theoretical calculations is an important tool for fullerene identification. However, the accuracy of available theoretical methods is often not adequate. Therefore, in this work, diﬀerent computational aspects needed to simulate realistically chemical shifts in the C 70 molecule are investigated by density functional theory (DFT) calculations. The importance of the functional choice, basis set, solvent, and molecular motions was assessed. The solvent was simulated using the implicit conductor-like polarized continuum model. The molecular motions were included via anharmonic corrections and averaging of snapshots obtained from classical and first-principles molecular dynamics (MD) simulations. Comparison to experiment revealed that density functional calculations typically overestimate the 13 C NMR chemical shifts. Hybrid functionals, such as BHandH and BHandHLYP, and long-range corrected functionals, such as wB97xd and CAM-B3LYP, give the best results. While the solvent has a minor eﬀect (chemical shift changes by B 1 ppm), the vibrational and dynamical eﬀects are surprisingly large, causing changes up to 9 ppm. Consideration of the latter was also necessary to explain the observed temperature dependence. While the dynamical corrections for MD performed in vacuum were overestimated, inclusion of the solvent in simulations provided more realistic results. The study thus points out the importance of an appropriate solvent model and a complex approach to the modelling, balancing the static, dynamic and environmental factors.


Introduction
Since the initial discovery of C 60 by Kroto et al. 1 the fullerene science has evolved into a rich multidisciplinary field.Fullerenes and their derivatives have been intensively studied as materials for solar cells, quantum computing, memory devices, magnetic storage devices, molecular containers, MRI contrast agents, and other applications. 2,3The fullerene cage structures have been rationalized by the isolated-pentagon rule (IPR). 4The number of the possible isomers of C n grows with n.For example, the C 60 fullerene has one isomer, C 80 has seven isomers, and for C 84 , 24 IPR isomers have been predicted to exist.In spite of the similarity of different fullerenes and their isomers, owing to their high molecular symmetry, they can be in many cases distinguished according to the number and intensity of spectral lines in 13 C nuclear magnetic resonance (NMR) spectra.The NMR spectroscopy, hand in hand with theoretical interpretation of the experimental NMR spectra, is thus an important analytical tool in fullerene chemistry.
][7][8][9][10][11] In a review by Heine, 12 however, typical accuracy of 13 C NMR DFT calculations was estimated to be 2 ppm, which may not be sufficient for larger fullerenes with a denser spacing (o1 ppm) of spectral lines.
To assess the reliability of commonly used methods for calculations of the 13 C chemical shift we use the C 70 fullerene as a convenient model.It contains five different symmetry-equivalent carbon types (assigned as C a -C e in Fig. 1) which cover the range of B130-150 ppm spanning to a large part typical fullerene 13 C chemical shifts. 13The performance of various density functionals and basis set convergence are analyzed.
In the next stage, we assess the effects of dynamics, temperature, and solvent on the calculated 13 C NMR chemical shifts.The effects of solvent are approximated using the conductor-like polarized continuum model (CPCM). 14Several approaches of involving the molecular motion in the computations are explored: averaging snapshots obtained from classical molecular dynamics (MD) and the first-principles molecular dynamics (FPMD) trajectories, and a quantum averaging of anharmonic vibrational corrections. 15,16
The basis set convergence was investigated using the IGLO-II, 40 IGLO-III, 40 and IGLO-IV 41,42 basis sets optimized for calculations of NMR parameters.A set of Jensen's polarization-consistent basis sets, pcS-n, was also tested. 43lculated isotropic nuclear shielding (s) was converted to chemical shift (d) using tetramethylsilane (TMS) as a reference and C 60 as a secondary reference, according to formula where s i is a signal of C 70 carbon and d(C 60,TMS ) = 143.15ppm. 44

Quantum vibrational averaging
The nuclear potential of C 70 was expanded to a Taylor series up to fourth powers of all normal mode coordinates Q i . 45Only cubic and semi-diagonal normal mode quartic constants (with two or more identical indices) were considered, as these can be obtained by a single normal mode numerical differentiation of harmonic force fields. 46The anharmonic force field was calculated at the Hartree-Fock (HF) level using the 6-31G basis set, while the shielding derivatives were obtained at the BP86/IGLO-II level.Gaussian 09 18 was used for the Hessian computations for geometries displaced in normal modes.Program S4 47 interfaced to Gaussian was used for the anharmonic vibrational averaging.
One of the limitations of anharmonic vibrational corrections performed here is the necessity to use Hartree-Fock approximation which overestimates harmonic frequencies by about 10%. 48,49hile DFT provides typically better harmonic frequencies deviating from experiment by only a few percents, [50][51][52][53] it was not affordable for such extensive calculations on C 70 .We do not suppose that the relatively small frequency error of the HF method would significantly change the anharmonic corrections added in a perturbation way only.
For vibrational ground state c n isotropic nuclear magnetic shielding was calculated as The s 1 and s 2 are the first and second normal mode isotropic shielding derivatives obtained numerically with a differentiation step of 0.2 Å.This parameter was found to be optimal for many systems including fullerenes before. 54Within the second-order degeneracy-corrected perturbational (PT2) scheme 46 the wave function is expanded in the harmonic oscillator basis. 16lternatively, the vibrational configuration interaction (VCI) was used for vibrationally averaged shieldings, where the Hamiltonian was diagonalized in the harmonic oscillator basis functions (|li, up to doubly excited states) limited by an interaction parameter |hl|V|ni/(E n À E l )| r 0.1 where n is a fundamental vibration or the ground state.
A third approximation tried was based on the harmonic oscillator limit for energies 55 and provided temperature-corrected shielding (TCS) as Although the anharmonic methods should be in principle equivalent, they differ because of the approximations used and the vast number of harmonic oscillator (HO) states that should be considered for C 70 .The PT2 method provides correction for a large number of HO states but may fail for large anharmonic effects and Fermi resonances.The resonances are better represented in VCI, where, however, a smaller number of the states can be considered.Neither the VCI nor PT2 implementation involves the temperature effects involved in TCS, where, however, simplified anharmonic energy and property surfaces are supposed.Therefore, we used all three methods to get a feeling of the significance of the anharmonic vibrational averaging of the chemical shift.

Molecular dynamics
The BP86/def-TZVP optimized structure of C 70 was minimized using the Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) procedure in Tinker.The minimized geometry was evolved for 1 ns.The classical MD trajectories with the MM3 force field 56 were obtained using the modified Beeman algorithm with a 1 fs integration time step at constant temperature T = 300 K.The MM3 force field was modified for the C 70 carbon atoms as described in the ESI.† Chemical shifts were calculated for 40 snapshots at the BHandHLYP/IGLO-III level and averaged.To minimize systematic error, the dynamical correction was calculated as a difference between the average and the shielding calculated for MM3-minimized geometry.

First-principles molecular dynamics
FPMD simulations were performed on the BP86/def-SVP Born-Oppenheimer potential energy surface.The leapfrog Verlet algorithm as implemented in the Turbomole code 17 was used to update the coordinates.To avoid problems due to the limited numerical precision in the density optimization and the force calculation, the thresholds for energy and CPHF equations convergence were tightened to 10 À8 a.u.A denser numerical grid denoted as ''m5'' was used.The time step was 80 a.u.
(1.93 fs).To achieve a better sampling of the configuration space we ran several microcanonical (NVE ensemble) dynamics trajectories 1.2 ps long.The initial velocities were randomly distributed and corresponded to the temperature of 300 K.A total of four MD trajectories of 1.2 ps were produced for the averaging of snapshots.The initial 0.2 ps part of each trajectory was discarded as an equilibration phase, and then 16 equidistant snapshots of each trajectory (64 in total) were taken for the 13 C NMR calculations (BHandHLYP/IGLO-III) and averaging.
To verify the time-correlation of the data the block-averaging method was used. 57The mean errors of average were found in the interval of AE0.2 to AE0.25 ppm for the NMR shielding of different carbons, which after subtraction from the C 60 reference give a statistical error of AE0.35 ppm in the dynamical corrections for shielding.The dynamical corrections were obtained as differences between the average and result for the BP86/def-SVP optimized structure.

NMR experiment
Fullerenes (C 60 and C 70 ) were purchased from Sigma. 13 C NMR spectra of C 70 were recorded in deuterated ortho-dichlorobenzene (ODCB) using a Bruker Advance IIt 600 and 500 MHz spectrometer equipped with a cryoprobe.The spectra were measured at several temperatures (278, 298, 318, 338, 358 and 378 K).A small amount of C 60 was detected (at 142.86 ppm) in the sample of C 70 .The chemical shifts were referenced to tetramethylsilane (TMS).

Basis set convergence
We used optimized BP86/def-TZVP C 70 structure 58 to investigate the effect of the basis set and DFT functional on the 13 C NMR chemical shifts in C 70 calculated for the molecule at rest in vacuo.The BP86/def-TZVP level previously provided results in an excellent agreement with the experimental structural parameters for the C 60 molecule. 58,59able 1 lists the 13 C chemical shifts for C 70 calculated using the different IGLO-n type basis sets.The IGLO-III results appear converged within 0.2 ppm as compared to the IGLO-IV basis set, which can be considered to be near the basis set limit.Corrections beyond IGLO-IV can be expected to lie below 0.2 ppm. 60Hence the error of the IGLO-III basis set calculations with respect to the basis set limit can be estimated with safe margins to be within 0.4 ppm.With a conservative error estimate of 0.4 ppm for 13 C shifts the IGLO-III basis set provides a good compromise between the accuracy and the computational cost.IGLO-IV is recommended if higher accuracy is needed.
Following the suggestion of a referee, the basis set convergence was also tested on the set of Jensen's polarization consistent pcS-n basis sets 43 known to provide accurate nuclear shieldings. 61he convergence of 13 C nuclear shieldings in C 60 obtained with pcS-n (n = 0, 1, 2, 3) can be seen in Fig. S1 (ESI †).Table S1 (ESI †) lists the 13 C chemical shifts for C 70 calculated using different pcS-n (n = 0, 1, 2) basis sets.It is apparent that the IGLO-n basis sets provided faster convergence than pcS-n.

Performance of various density functionals
The results obtained using the selected DFT functionals are summarized in Table 2, where the calculated 13 C shifts as well as their root-mean-square deviations (RMSD) from the experimental values are collected.The results in Table 2 are ordered according to the RMSD.Typically, the DFT calculations overestimate the experimental 13 C NMR chemical shift by several ppm.Obviously, the effects of solvent and temperature are not included at this stage, and we can compare only the main trends exhibited by the functionals.
Nevertheless, for the gas-phase C 70 molecule at rest it appears that the long-range corrected functionals CAM-B3LYP and wB97xd as well as the hybrid functionals with 50% of exactexchange admixture BHandH and BHandHLYP perform the best, while the non-hybrid GGA functionals give a worse agreement with the experiment.The OPBE functional, provides results of the same accuracy as BHandHLYP with much smaller CPU cost.5][66][67] A good performance is provided also by the M06L pure DFT functional of Truhlar and Zhao. 33ombining different exchange functionals with the LYP correlation along the B1LYP, B3LYP, X3LYP, and O3LYP series shows rather marginal dependence of the performance on the exchange functional.Increasing the exact exchange admixture in the functional generally lowers the calculated 13 C chemical shifts, bringing them closer to the experimental ones, as seen in Fig. 2. A decrease of the RMSD with increasing exact-exchange admixture is observed for relative carbon shifts (C signals relative to C e ) in Fig. 2. Noticeably, if one looks at relative 13 C NMR chemical shifts in Table S2 (ESI †) instead of the absolute ones in Table 2 the relative performance is somewhat altered especially for larger RMDS values.
0][71][72][73][74][75] However, less significant improvement for carbon only containing compounds was observed. 76f the larger IGLO-IV basis set was used in the calculations shown in Table 2 (not affordable presently), most RMSDs should improve, as the DFT calculated 13 C chemical shifts would lower by about 0.1-0.2ppm, mostly towards the experimental values.We do not expect this rather small basis set effect to alter the performance of the DFT functionals, as shown in Table 2.

The solvent effect
As seen in Table 3, the inclusion of the solvent within the CPCM implicit solvent model moderately decreases the calculated 13 C shifts by 0.5-1.0ppm towards the experimental values.Note that the new experiment performed in this work in ODCB solvent (e = 9.93) Fig. 2 The effect of increasing exact-exchange admixture in DFT functionals on calculated 13 C chemical shifts.RMSD (solid line) stands for root mean square deviations from the experimental values. 13RMSD r (dashed line) is calculated for relative 13 C shifts (C-C e ).  1).This suggests a limited dependence of the 13 C NMR chemical shifts on the kind of solvent for this non-polar molecule.A closer look at Table 3 reveals that the BHandH and OLYP functionals decrease calculated shifts towards experimental results (RMSD drop of B0.1 ppm), while including the effects of solvent slightly increases the RMSD error at the wB97xd level, from 0.4 to 0.5.

Equilibrium geometry effect
The effect of the equilibrium geometry on the calculated 13 C shifts for C 70 is illustrated in Table 4.While improving the basis set from def-TZVP towards def2-QZVP, the shift decreases towards the experiment in a uniform way, by about 1 ppm, when comparing the def-TZVP and def2-QZVP structures.
Changing the DFT functional used for the optimization leads to substantially larger shift changes (obtained at the BHandHLYP/IGLO-III level), and they can both improve or worsen the results, as seen in Table 4.The B3LYP structure provides chemical shifts closer to experiment, while the B97d structure largely overestimates the shifts.If relative chemical shifts are compared to experimental results in Table 4 via RMSD r (obtained from C signals relative to C e ), the B3LYP functional provides worse agreement than BP86 and B97d.
The seemingly good performance of, e.g. the B3LYP/def2-TZVP structure as compared to our default level BP86/def-TZVP is presumably due to error cancellation, as seen in Fig. 3, where the C-C distances in C 60 calculated at different levels are compared to the experimental ones.The B3LYP/def2-TZVP or BHandHLYP/def2-TZVP levels provide worse calculated C-C distances than our default level for structure optimization BP86/def-TZVP.The BP86/def-TZVP level appears to provide reliable molecular structures also for endohedral helium fullerenes. 58,59e molecular motions and 13

C chemical shifts
The results summarizing different strategies for accounting of the temperature and dynamics are summarized in Table 5, as corrections to individual carbon nuclear magnetic shieldings, i.e. differences between dynamically/vibrationally averaged and optimized-structure nuclear shieldings.Averaged distances and spatial fluctuations of carbon atoms derived from the simulations are detailed in Table S3 (ESI †).
The changes caused by the molecular motions for the gasphase C 70 molecule can be assigned as deshielding and are relatively large, up to ca.À9 ppm, very much depending on the method.The positive (shielding) contribution result for C e in the MD line in Table 5 may be an artifact caused by the inaccurate force field.The vibrational averaging in PT2 approximation gives nearly uniform deshielding of about À4.5 ppm for all carbon types, while the TCS, classical MD, and FPMD approaches show relatively more variable effects.The corrections for carbon C e are significantly different from the other carbon types.The geometric fluctuations of atoms C a -C e around the equilibrium position are rather uniform (Table S4, ESI †), hence we infer that the differential dynamical contributions for C a -C e in FPMD originate from the fluctuations in the electronic structure.
For the C 60 reference, dynamical correction to 13 C nuclear magnetic shielding about À3 to À2 ppm was obtained using different methods comparable with that of À2.5 ppm mentioned previously. 78If anomalous contribution for atom C e in C 70 is    2 are added to the dynamical corrections in Table 5 (taken with opposite sign to turn the nuclear magnetic shielding into the chemical shift), the agreement with the experiment is not improved.However, when the FPMD simulations are performed in the solvent, the dynamical corrections (Table 5) for the 13 C nuclear magnetic shielding are still deshielding but now considerably smaller in absolute size and more uniform (À1.2 to À1.7 ppm) than those obtained for FPMD simulation in vacuo.A detailed look reveals that the solvent damps the molecular motions.This is documented in Table S2 (ESI †) on the average deviations from equilibrium positions.
The different contributions to the total 13 C NMR chemical shifts in C 70 discussed above are summarized in Fig. 4. The inclusion of solvent and basis set correction decreases the gas phase shifts and thus improve the agreement with experiment.
The contributions due to molecular motions (FPMD is shown as example in Fig. 4) are mostly positive, and shift the values away from the experimental ones.If the FPMD simulations could be run longer, an improvement of a couple of tenths of ppm could be expected due to somewhat large statistical error (mean error of the average, see the Methods section) obtained, AE0.35 ppm for the dynamical corrections to the 13 C chemical shift.The total calculated 13 C shifts corrected for the basis set error and dynamical solvent are again closer to the experimental ones.However, the inclusion of the dynamical corrections including solvent does not bring any significant improvement over the static molecule in vacuo results, because they are small in size, relatively uniform, and largely cancelled by the dynamical effect of the reference C 60 compound.

The temperature dependence
Vibrational averaging obtained using the TCS method also enables the estimation of the temperature dependence of the carbon chemical shifts in C 70 .Results for 273 K and 373 K calculated using two different computational levels are shown in Table S4 (ESI †).The calculated and experimental changes (related to value for carbon C e ) are shown in Fig. 5.The temperature dependence appears to be relatively uniform for all carbons at both levels, B0.6-0.9 ppm per 100 K. Experimental dependence can be seen in Fig. S2 (ESI †).The correlation between experimental and theoretical results is clear.However, the observed change is about 5 times smaller than in theory.This can be partially attributed to the missing solvent effects in the TCS calculations, which, however, could not be included due to the computational time limits.

Conclusions
To estimate the accuracy of available theoretical procedures, 13 C NMR chemical shifts in the C 70 fullerene were modelled using various density functionals and basis sets.The effects of solvent, dynamics, and temperature on the 13 C shifts in the theoretical simulations were discussed.Reasonable 13 C chemical shift values within several ppm from the experimental results could be obtained already by the usual calculations of a gasphase molecule at rest.With hybrid or long-range corrected functionals the calculated chemical shifts were overestimated by about 2 ppm.The largest errors, up to 5 ppm, were found for the GGA functionals, such as BP86.The basis set effects on the 13 C shifts calculated for the IGLO-II to IGLO-IV series are moderate; the errors using the IGLO-III basis set are estimated within 0.4 ppm of the basis set limit, or 0.2 ppm from IGLO-IV results.Improving the basis set size in the structure optimization resulted in changes a few tenths of ppm on the 13 C NMR  chemical shifts, whereas changing the functional had a larger impact of several ppm.
Inclusion of the implicit solvent effects caused a decrease of the calculated 13 C chemical shifts towards the experimental values by about 1 ppm.The effect of molecular motions on the 13 C NMR chemical shift in gas-phase C 70 caused relatively large deshielding contributions of about 2-9 ppm, depending on the method used for averaging.These changes were strongly damped if the solvent was considered in simulation, leading to a deshielding contribution of about 1.2-1.7 ppm only.For the chemical shifts, the changes were largely cancelled by the motional contributions from the C 60 reference.

Fig. 3
Fig. 3 Calculated bond lengths in optimized C 60 obtained using the BP86 functional and different basis sets (on the right) and density functionals in combination with the def2-TZVP basis set (on the left).Experimental results are obtained from ref. 77.The C-C bonds are displayed on the top.

Fig. 4
Fig. 4 Importance of various effects on calculated 13 C chemical shifts in C 70 .BHandHLYP stands for shifts obtained by the vacuum BHandHLYP/IGLO-III calculation, ''solvent'' denotes the BHandHLYP/IGLO-III/CPCM model, and the other bars are related to the involvement of vacuum and solvent FPMD dynamics.All data are corrected for the estimated basis set error calculated from difference between the IGLO-IV and IGLO-III results in Table1(all results comprise the previous enhancements, e.g. the ''dynamic'' = BHandHLYP shifts in vacuum + solvent contribution + basis set correction + dynamic contribution).The experiment is taken from ref.13.

Fig. 5
Fig. 5 Temperature effect on NMR shifts in C 70 .Plot shows comparison of experimental and theoretical (BP86/IGLO-II//HF/6-31G) differences in relative shieldings caused by temperature (D = 100 K).The shift is referenced to C e .

Table 1
13sis set convergence of13C chemical shifts (in ppm) for C 70

Table 3
Vacuum vs. solvent chemical shifts (in ppm) calculated in C 70 provided almost identical results to the previous experiments done in 1,1,2,2-tetrachloroethane (see Table

Table 4
13lculated (BHandH/IGLO-III)13C chemical shift (in ppm) for C 70 using geometries optimized at different computational levels The def2-TZVP basis set was used.RMSD is the root-mean-square deviation from experimental values ref.13, RMSD r is calculated for relative chemical shifts (C-C e ).
a The BP86 functional was used.b

Table 5
13namical effects on13C nuclear magnetic shieldings (ppm) presented as difference between vibrationally averaged and minimum structure results excluded then the ''dynamical'' contribution for C 60 represents approximately 70% of average C 70 values (MD, FPMD).The vibrational (PT2, VCI, TCS) contributions to the13C NMR chemical shifts in C 70 are likely to be similar as in C 60 .Therefore, the usage of C 60 as the secondary reference will partly cancel the dynamical corrections.If the best DFT results in Table