October 01, 2024
Theoretical Study of Reactivity Indices and Rough Potential Energy Curves for the Dissociation of 59 Fullerendiols in Gas-Phase and in Aqueous Solution with an Implicit Solvent Model
Anne Justine ETINDELE *Higher Teachers Training College, University of Yaounde I, P.O. Box 47, Yaounde, CAMEROON
e-mail: annetindele@yahoo.fr*
Abraham PONRA
*Department of Physics, Faculty of Science, University of Maroua, P.O.Box 814, Maroua, CAMEROON
e-mail: abraponra@yahoo.com*
Mark E. CASIDA
*Laboratoire de Spectrométrie, Interactions et Chimie théorique (SITh), Département de Chimie Moléculaire (DCM, UMR CNRS/UGA 5250), Institut de Chimie Moléculaire de Grenoble (ICMG, FR2607), Université Grenoble Alpes (UGA) 301 rue de la Chimie, BP 53,
F-38041 Grenoble Cedex 9, FRANCE
e-mail: mark.casida@univ-grenoble-alpes.fr*
G.Andrés CISNEROS
*Department of Physics, Department of Chemistry and Biochemistry, University of Texas at Dallas, Richardson, Texas, 75802, USA
e-mail: andres@utdallas.edu*
Jorge NOCHEBUENA
*Department of Chemistry and Biochemistry, College of Science and Mathematics, Augusta University, Augusta, Georgia, 30912, USA
e-mail: jnochebuenaherna@augusta.edu*
Abstract
Buckminsterfullerene, C\(_{60}\), has not only a beautiful truncated icosahedral (soccerball) shape, but simple Hückel calculations predict a three-fold degenerate lowest unoccupied molecular orbital (LUMO) which can accomodate up to six electrons making it a good electron acceptor. Experiments have confirmed that C\(_{60}\) is a radical sponge and it is now sold for use in topical cosmetics. Further medical uses require functionalization of C\(_{60}\) to make it soluble and one of the simplest functionalization is to make C\(_{60}\)(OH)\(_n\) fullerenols. A previous article [Adv.Quant.Chem. 8, 351 (2023)] studied reactivity indices for the successive addition of the \(^\bullet\)OH radical to (\(^\bullet\))C\(_{60}\)(OH)\(_n\) in gas phase. [(\(^\bullet\))C\(_{60}\)(OH)\(_n\) is only a radical when \(n\) is an odd number.] This present article extends this previous work by examining various aspects of how the reaction, \[^\bulletC_{60}OH + ^\bulletOH \rightarrow C_{60}(OH)_2 \, \, \, (1) \nonumber\] changes in aqueous solution. One obvious difference between C\(_{60}\) and their various isomers of C\(_{60}\)(OH)\(_2\) is the presence of a dipole. As fullerendiols are nearly spherical, their change in dipole moment in going from gas to aqueous phase may be estimated using back-of-the-envellope calculations with the Onsager model. The result is remarkably similar to what is obtained using density-functional theory (DFT) witn an implicit solvation model (Surface Molecular Density, SMD). Calculation of fullerendiol C-O bond energies and reactivity indices using the SMD approach confirm that the general conclusions from the earlier work regarding gas-phase reactivity still hold in the aqueous phase. A major difference between the present work and the earlier work is the calculation of potential energy curves (PECs) for reaction (1) in gas and aqueous phases. This is done in exploratory work for all 59 possible fullerendiols in both gas phase and in aqueous solution with the SMD approach using spin-unrestricted DFT calculations with symmetry breaking. Surprisingly little change is found between the gas-phase and aqueous-phase PECs. It is discovered that the majority of C\(_{60}\)(OH)\(_2\) show radicaloid character, as might have been expected from trying to draw resonance structures. Spin-contamination curves are also remarkably similar for gas-phase and aqueous-phase results. Although our calculations do not include a dispersion correction, it was noticed that all calculated PECs have a \(1/R^6\) behavior over a significant \(R = R(C-O)\) distance, underlying the need to be careful of double counting when including dispersion corrections in DFT. A short coming of our SMD approach is the lack of explicit water molecules which can form hydrogen bonds with the OH groups and dissociating radicals.
Buckminsterfullerene C\(_{60}\) is well-known for its geometrical beauty, reminiscent of the geodesic domes of Buckminsterfuller. It can also capture up to six electrons in its \(t_{1u}\) lowest unoccupied molecular orbitals (LUMOs). This strongly electrophilic property has been characterized experimentally and it has been called a “radical sponge” [1]. It is even used as a commercial ingredient in some skin care products, one of which is an anti-aging moisturizer which goes by the name “C60” [2]. However many applications require an increase in the aqueous solubility of this hydrophobic molecule. A particularly simple way to increase the solubility of C\(_{60}\) is to decorate it with hydroxyl groups. As C\(_{60}\) is a “radical sponge,” it can react with a large number of hydroxyl radicals [3]–[16] to create fullerenols (\(^\bullet\))C\(_{60}\)(OH)\(_n\) where the bullet (\(^\bullet\)) is a reminder that these fullerenerols are radicals for odd \(n\).
To take into account the impact of solvents on the properties of molecules, two approaches may be considered. The first approach takes the molecular nature of the solvent molecules explicitly into account. We will only use the second in the present work. This is the implicit model in which the solvent is considered as a dielectric continuum. As the (\(^\bullet\))C\(_{60}\)(OH)\(_n\) are nearly spherical, Onsager’s spherical dielectric cavity reaction field model provides a pencil and paper way to determine the dipole moment for the molecule in solution from its gas-phase dipole moment. Thus Onsager’s spherical dielectric cavity reaction field model provides a pencil and paper way to determine the dipole moment for the molecule in solution from its gas-phase dipole moment. This first-order approximation works almost surprisingly well compared to other implicit solvent models that use cavities which are more carefully adjusted to reflect the shape of the molecule.
The global and local reactivity indices can be significantly impacted depending on whether the study is carried out in the gas or aqueous phase [17], [18]. These reactivity indices made it possible to effectively study the donor and acceptor character of C\(_{60}\)(OH)\(_n\) [16] in the gas phase. So, what about this character for the aqueous phase? We are given the opportunity to try to answer this question in this paper. In addition, the description of potential energy curves (PECs) could provide valuable information on the impact of the presence of the solvent on the dissociation (or formation) of the fullenerol-hydroxyl (i.e., \(^\bullet\)C\(_{60}\)(OH) + \(^\bullet\)OH) bonds.
Previous work has shown that the dissociation of a molecule into two fragments requires taking into account symmetry breaking along the PECs [19], [20]. Taking into account this symmetry breaking, point-by-point, throughout the PECs provides a unique opportunity to determine the Coulson-Fisher points for the 59 isomers of fullerenediol by examining the variation of \(\langle {\hat{S}}^2 \rangle\) with distance.
This article is organized in the following manner: The next section covers the basic theory used in the rest of the paper. Section 3 gives the computational details that we are using in this paper for all the calculations. Section 4 presents and discusses our results. In particular, we present PECs and their associated \(\langle {\hat{S}}^2 \rangle\) for all 59 isomers in the SI. Section 5 is the concluding discussion.
Electronic structure calculations were initially done for atoms, gradually moved onto small molecules, and then—with the development of band theory—onto crystals. However it is important to be able to model the electronic structure of solvated molecules for the simple reason that most experiments are carried out in solution, or at least involve key steps and/or measurements carried out in solution. We use an implicit solvent model. This replaces the solvent with a dielectric cavity. It has the advantage of being relatively simple and so computationally expedient, but does not take into account all of the physics of a molecule in solution.
The implicit solvent model used in the present work is the well-known SMD (for Solvation Model based upon the quantum mechanical Density) [21]. This is a reaction field model where the solvent molecules are replaced by a dielectric continuum surrounding the solvated molecules in a solvent excluded cavity (SEC, also said to be defined by a solvent accessible surface). The molecule in the cavity is described by quantum mechanics (QM), from which its charge distribution is calculated. This charge distribution creates a reaction field in the surrounding dielectric which, in turn, modifies the QM calculation. Although this sounds like it requires an iterative self-consistent calculation, its solution may be done in a non-iterative single-step procedure. Particular implicit solvent models vary by how the SEC is generated and how the Poisson equation is solved to determine the reaction field. In particular, the generation of the SEC involves some semi-empirical parameters which are fit, in the SMD model, to experimental solvation free energies. It is not our purpose to go further into details here.
However fullerenediols are nearly spherical molecules which means that the oldest most approximate reaction field model—namely the Onsager model [22] is expected to be a good approximation. Although a program such as Gaussian [23] has an option to do Onsager reaction field calculations, it is really a shame to use such a large program for what is ultimately a back-of-the-envelope calculation. Indeed, doing the back-of-the-envelope calculation leads to enhanced understanding and so we do present this model here. As with more sophisticated implicit solvent models, the solvent is treated as a polarizable continuum with dielectric constant \(\epsilon\). However the molecule is treated as a point dipole moment \(\mu_{\ell}\) whose value will depend upon its surroundings. It also has a polarizability \({\boldsymbol{\alpha}}\). The molecule is placed in a solvent-excluded spherical cavity of radius \(a\). The electric field of the molecular dipole moment will cause a redistribution of charges around the sphere and induce a reaction (electric field) \({\vec{\cal E}}\) inside the sphere which acts on the molecule to change the dipole moment through the equation, \[{\vec{\mu}}_{\ell} = {\vec{\mu}}_g + {\boldsymbol{\alpha}} {\vec{\cal E}} \, , \label{eq:Onsager461}\tag{1}\] where \({\vec{\mu}}_{\ell}\) is the molecular dipole moment in the liquid and \(\mu_g\) is the gas-phase dipole moment. This reaction field is (Ref. [24] pp. 130-134) \[{\vec{\cal E}} = \frac{2(\epsilon-1)}{2 \epsilon+1} \frac{1}{a^3} {\vec{\mu}}_{\ell} \, . \label{eq:Onsager462}\tag{2}\] Therefore, \[{\vec{\cal E}} = \frac{2(\epsilon-1)}{2 \epsilon+1} \frac{1}{a^3} \left( {\vec{\mu}}_g + {\boldsymbol{\alpha}} {\vec{\cal E}} \right) \, , \label{eq:Onsager463}\tag{3}\] which can be solved for the reaction field to give, \[{\vec{\cal E}} = \left( \frac{2\epsilon+1}{2(\epsilon-1)} a^2 {\boldsymbol{1}} - {\boldsymbol{\alpha}} \right)^{-1} {\vec{\mu}}_g \, . \label{eq:Onsager464}\tag{4}\] Finally, plugging this result into Eq. (1 ) and rearranging gives, \[{\vec{\mu}}_{\ell} = \left( {\boldsymbol{1}} - \frac{2(\epsilon-1)}{2\epsilon+1} \frac{1}{a^3} {\boldsymbol{\alpha}} \right)^{-1} {\vec{\mu}_g} \, . \label{eq:Onsager465}\tag{5}\] This still leaves the problem of the polarizability and the size of the cavity. However assuming an isotropic medium, so that we may replace the polarizability tensor \({\boldsymbol{\alpha}}\) with the average polarizability \({\bar \alpha}\), gives, \[{\vec{\mu}}_{\ell} = \left( 1 - \frac{2(\epsilon-1)}{2\epsilon+1} \frac{\bar \alpha}{a^3} \right)^{-1} {\vec{\mu}_g} \, , \label{eq:Onsager466}\tag{6}\] so that we actually only need to determine \({\bar \alpha}/a^3\). This can be done using the Lorenz-Lorentz equation relating the polarizabilty and the refractive index \(n\), \[\frac{n^2-1}{n^2+2} = \frac{4\pi}{3} \rho {\bar \alpha} \, , \label{eq:Onsager467}\tag{7}\] where \(\rho\) is the number density so that, \[\frac{1}{\rho} = \frac{4}{3} \pi a^3 \, . \label{eq:Onsager468}\tag{8}\] Hence, \[\frac{\bar \alpha}{a^3} = \frac{n^2-1}{n^2+2} \, . \label{eq:Onsager469}\tag{9}\] Note that using these equations for treating C\(_{60}\) in water requires the dielectric constant for water but the refractive index of C\(_{60}\). For water, \(\epsilon=80\) at 20 \(^\circ\)C (Ref. [25], p. 299). The refractive index of solid C\(_{60}\) which is 2.2 at 630 nm wavelength [26]. Putting everything together, then \[{\vec{\mu}}_{\ell} = 2.23 {\vec{\mu}}_g \, . \label{eq:Onsager4610}\tag{10}\] for C\(_{60}\) in aqueous solution. We expect similar results for fullerenols.
The International Union of Pure and Applied Chemistry (IUPAC) numbering (Fig. 1) is used throughout for C\(_{60}\) [27]. An “xyz” format file with this numbering has been supplied in the Electronic Supplementary Information (ESI) for use with molecular visualization software.

Figure 1: IUPAC numbering of buckminsterfullerene [27]..
The same level of electronic structure calculation is used throughout this article, whether it be for the gas-phase or for the aqueous phase using an implicit solvent model.
Calculations were carried out using version 5.0.4 of Orca [28] and with Gaussian [23]. The B3LYP keyword was used with Gaussian while the B3LYP/g keyword was used with Orca. This is the same functional which we will rename B3LYP(VWN3) [29], [30] because the original B3LYP functional (programmed in Gaussian) used the Vosko-Wilk-Nusair parameterization of the random-phase approximation correlation energy (known as VWN3 in Gaussian) but that the default B3LYP functional in Orca uses the Vosko-Wilk-Nusair parameterization of the Ceperley and Alder’s quantum Monte Carlo correlation energy (known as VWN5 in Gaussian) unless “/g” is specified. The same def2SVP [31] orbital basis set was used in all our calculations.
The default in Orca is also to use a resolution-of-the-identity (RI) approximation to reduce the problem of calculating four-center electron repulsion integrals to that of calculating only three-center electron repulsion integrals. As this option was used in our previous gas-phase work [16] on gas-phase reactivity indices, we continued with the same option in the article when calculating reactivity indices in solution. It was also used in calculating gas and solution phase dipole moments. However it was judged useful to follow the advice of Ref. [32] in order to ensure the best possible agreement between Orca and Gaussian calculations:
” In ORCA, keywords used included
NORI(no approximation is used),TightSCF,TightOpt,SlowConv,NumFreqand an ultra-fine grid. We explicitly confirmed that Gaussian and Orca gave the same results for the same basis sets and functionals in single-point calculations.” [32]
Implicit solvent calculations were carried out with Orca [28] at the B3LYP(VWN3)/def2-SVP level using the SMD [21] solvent model. It has been discovered, in work on polypyridine ruthenium clusters, that Gaussian and Orca use different grids on the solvent-allowed surface which leads to significant differences in numerical results (see the ESI of Ref. [33]). Following that work, we opted to use
num_leb 590
end
in the Orca input to in order ensure good agreement with Gaussian.
| Count | Symmetry-Equivalent Isomers | Distance |
|---|---|---|
| 2 | (1,2)-C\(_{60}\)(OH)\(_2\), (1,5)-C\(_{60}\)(OH)\(_2\) | 1 bond |
| +2 = 4 | (1,3)-C\(_{60}\)(OH)\(_2\), (1,4)-C\(_{60}\)(OH)\(_2\) | 2 bonds |
| +2 = 6 | (1,6)-C\(_{60}\)(OH)\(_2\), (1,12)-C\(_{60}\)(OH)\(_2\) | 2 bonds |
| +2 = 8 | (1,7)\(^*\)-C\(_{60}\)(OH)\(_2\), (1,11)\(^*\)-C\(_{60}\)(OH)\(_2\) | 3 bonds |
| +2 = 10 | (1,8)\(^*\)-C\(_{60}\)(OH)\(_2\), (1,10)\(^*\)-C\(_{60}\)(OH)\(_2\) | 2 bonds |
| +1 = 11 | (1,9)-C\(_{60}\)(OH)\(_2\) | 1 bond |
| +2 = 13 | (1,13)-C\(_{60}\)(OH)\(_2\), (1,20)-C\(_{60}\)(OH)\(_2\) | 3 bonds |
| +2 = 15 | (1,14)-C\(_{60}\)(OH)\(_2\), (1,19)-C\(_{60}\)(OH)\(_2\) | 4 bonds |
| +2 = 17 | (1,15)\(^*\)-C\(_{60}\)(OH)\(_2\), (1,18)\(^*\)-C\(_{60}\)(OH)\(_2\) | 3 bonds |
| +2 = 19 | (1,16)-C\(_{60}\)(OH)\(_2\), (1,17)-C\(_{60}\)(OH)\(_2\) | 4 bonds |
| +2 = 21 | (1,21)-C\(_{60}\)(OH)\(_2\), (1,30)-C\(_{60}\)(OH)\(_2\) | 4 bonds |
| +2 = 23 | (1,22)\(^*\)-C\(_{60}\)(OH)\(_2\), (1,29)\(^*\)-C\(_{60}\)(OH)\(_2\) | 4 bonds |
| +2 = 25 | (1,23)-C\(_{60}\)(OH)\(_2\), (1,28)-C\(_{60}\)(OH)\(_2\) | 5 bonds |
| +2 = 27 | (1,24)\(^*\)-C\(_{60}\)(OH)\(_2\), (1,27)\(^*\)-C\(_{60}\)(OH)\(_2\) | 4 bonds |
| +2 = 29 | (1,25)-C\(_{60}\)(OH)\(_2\), (1,26)-C\(_{60}\)(OH)\(_2\) | 3 bonds |
| +2 = 31 | (1,31)-C\(_{60}\)(OH)\(_2\), (1,40)-C\(_{60}\)(OH)\(_2\) | 5 bonds |
| +2 = 33 | (1,32)\(^*\)-C\(_{60}\)(OH)\(_2\), (1,39)\(^*\)-C\(_{60}\)(OH)\(_2\) | 6 bonds |
| +2 = 35 | (1,33)\(^*\)-C\(_{60}\)(OH)\(_2\), (1,38)\(^*\)-C\(_{60}\)(OH)\(_2\) | 5 bonds |
| +2 = 37 | (1,34)\(^*\)-C\(_{60}\)(OH)\(_2\), (1,37)\(^*\)-C\(_{60}\)(OH)\(_2\) | 5 bonds |
| +2 = 39 | (1,35)-C\(_{60}\)(OH)\(_2\), (1,36)-C\(_{60}\)(OH)\(_2\) | 6 bonds |
| +2 = 41 | (1,41)-C\(_{60}\)(OH)\(_2\), (1,48)-C\(_{60}\)(OH)\(_2\) | 6 bonds |
| +2 = 43 | (1,42)-C\(_{60}\)(OH)\(_2\), (1,47)-C\(_{60}\)(OH)\(_2\) | 6 bonds |
| +2 = 45 | (1,43)-C\(_{60}\)(OH)\(_2\), (1,46)-C\(_{60}\)(OH)\(_2\) | 6 bonds |
| +2 = 47 | (1,44)-C\(_{60}\)(OH)\(_2\), (1,45)-C\(_{60}\)(OH)\(_2\) | 5 bonds |
| +2 = 49 | (1,49)-C\(_{60}\)(OH)\(_2\), (1,55)-C\(_{60}\)(OH)\(_2\) | 7 bonds |
| +2 = 51 | (1,50)-C\(_{60}\)(OH)\(_2\), (1,54)-C\(_{60}\)(OH)\(_2\) | 7 bonds |
| +2 = 53 | (1,51)-C\(_{60}\)(OH)\(_2\), (1,53)-C\(_{60}\)(OH)\(_2\) | 6 bonds |
| +1 = 54 | (1,52)-C\(_{60}\)(OH)\(_2\) | 8 bonds |
| +2 = 56 | (1,56)\(^*\)-C\(_{60}\)(OH)\(_2\), (1,59)\(^*\)-C\(_{60}\)(OH)\(_2\) | 8 bonds |
| +2 = 58 | (1,57)-C\(_{60}\)(OH)\(_2\), (1,58)-C\(_{60}\)(OH)\(_2\) | 7 bonds |
| +1 = 59 | (1,60)-C\(_{60}\)(OH)\(_2\) | 9 bonds |
A thorough study of the gas-phase reactivity of fullerenols with respect to step-wise addition of hydroxyl radicals has already been presented in Ref. [16]. Our interest here is in how the physical and chemical properties of fullerenols differ in gas phase and in aqueous solution. We will confine our study to fullerenediols only. As Table 1 shows, this is still a very large number of isomers. We will begin with physical properties and then chemical reactivity indices before going on to potential energy curves (PECs) for the dissociation reaction, \[C_{60}(OH)_2 \rightarrow ^\bulletC_{60}OH + ^\bulletOH \, . \label{eq:results461}\tag{11}\]

Figure 2: Graph correlating the magnitudes of gas-phase and SMD aqueous-phasedipole moments for C\(_{60}\)(OH)\(_2\) fullerenols calculated usinggeometries optimized in the gas phase. The dashedline is the result of Onsager’s implicit solvent model. See text..
SMD implicit solvent calculations of the dipole moments of the fullerenediols in solution were carried out at the spin-restricted (same-orbitals-for-different-spins or SODS) level using geometries optimized in the gas phase. The simplest implicit solvent model is the Onsager model which we expect to be a good first approximation for fullerendiols. This is confirmed in Fig. 2 where we present a correlation plot between calculated gas-phase and aqueous-phase dipole moments. The Onsager model predicts a straight line passing through the origin and with slope, \[m = 1 - \left( \frac{2 (\epsilon -1)}{2 \epsilon+1} \right) \left( \frac{n^2-1}{n^2+2} \right) \, . \label{eq:SMD461}\tag{12}\] Here \(\epsilon\) is the dielectric constant of the surrounding water, which is well-known to be close to 80, and \(n\) is the refractive index of solid C\(_{60}\) which is 2.2 at 630 nm wavelength [26]. Using these values gives a slope of \(m=0.44906\) which corresponds to the dashed line in Fig. 2. We can conclude that the Onsager model is a reasonable first approximation to the SMD implicit solvent model in this case.
Given such a large change in the dipole moments going between the gas- and aqueous-phases, we might also except a large change in reactivity indices. As in Ref. [16], we seek how well reactivity indices predict the C-O bond dissociation energy (BDE). This latter is calculated as the difference between the sum of the spin-unrestricted (different-orbital-for-different-spins or DODS) energies of \(^\bullet\)C\(_{60}\)(OH) + \(^\bullet\)OH and the SODS energy of C\(_{60}\)(OH)\(_2\). The gas-phase results confirmed that the BDE energy increases (i.e., the fullerenediol becomes more stabe) as the spin density of the reaction site \(^\bullet\)C\(_{60}\)(OH) increases. Part (a) of Fig. [fig:SMDReactivityIndices] shows that this also holds in our solution study. The gas-phase results also confirmed the electrophilic nature of the \(^\bullet\)OH radical because larger BDEs were generally found at sites with more negative Mulliken charges. Part (b) of Fig. [fig:SMDReactivityIndices] shows that this continues to hold in aqueous solution. Finally large values of the radical Fukui function \(f^0\) were also found to favor large values of the BDE in the gas phase and part (c) of Fig. [fig:SMDReactivityIndices] shows that this continues to hold in aqueous solution. That the SMD aqueous solution results look very much like the corresponding gas-phase results [16] is a nice confirmation that \(^\bullet\)OH remains an electrophilic radical in solution even with respect to an “electron sponge” such as C\(_{60}\).
The large number of fullerenediol isomers (not to mention their various conformers!) shown in Table 1 represents a serious complication for studying the PECs for reaction 11 . In order to keep things manageable, we have chosen to only calculate PECs for for rigid vertical removal of \(^\bullet\)OH without any internal relaxation of either \(^\bullet\)C\(_{60}\)OH or \(^\bullet\)OH (Fig. 4).

Figure 4: Illustration for isomer 1 of a procedure used for isomers 1,2, and 3. The vertical docking of \(^\bullet\)OH on\(^\bullet\)C\(_{60}\)OH was accomplished by begining from the fullerendioland using Molden [34] to construct a Z-matrix. Then onlythe C-O distance was varied for carbon number 1, keeping all other bondand dihedral angles fixed. In particular the structures of \(^\bullet\)OHand \(^\bullet\)C\(_{60}\)OH were kept fixed..
Calculations were of the spin-unrestricted (DODS) type with symmetry breaking so as to be able to describe the breaking of the C-O single bond as correctly as possible. Both energies and spin contamination (\(\langle \hat{S}^2 \rangle\)) were monitored. This produced much too many results to given in the main body of the paper so we will just give a few examples and some summarizing statistics. Full results may be found in the Supplementary Information (SI) associated with this article. What we expected to see were curves such as those found for 1,5-C\(_{60}\)(OH)\(_2\) and shown in Fig. 5. At large \(R\)(C-O), we expect to see a heavily spin-contaminated open-shell singlet with \(\langle \hat{S}^2 \rangle \approx 1\). As the \(^\bullet\)C\(_{60}\)(OH) and \(^\bullet\)OH radicals approach and form a bond, we expected a normal closed-shell singlet with \(\langle \hat{S}^2 \rangle = 0\). Interestingly, this only happened for about a quarter of the symmetry pairs studied.

Figure 5: 1,5-C\(_{60}\)(OH)\(_2\) (isomer 6) vertical dissociation PECsand spin-contamination for the same geometries in gas phase and inaqueous solution..
Much like with the reactivity indices, the PECs and spin-contamination curves turn out to be remarkably similar in gas-phase and in aqueous solution. The degree of spin contamination for the different isomers is shown in Fig. 7. It is apparent that there are many fullerediols which are not simple closed-shell singlets. At first thought, this may be a bit surprising as such structures are rarely mentioned in the relevant literature. On second thought, it is easy to understand why radicaloid structures might be preferred in many cases. One such structure is illustrated in Fig. 6 and has been verified by examination of the spin density of 1,8-C\(_{60}\)(OH)\(_2\). In contrast, both 1,7-C\(_{60}\)(OH)\(_2\) and 1,9-C\(_{60}\)(OH)\(_2\) are closed-shell singlets.

Figure 6: Full (top) and partial (bottom) fullerendiol structures showingthe diradical nature of 1,8-C\(_{60}\)(OH)\(_2\). Note that there areseveral other resonance structures that could be drawn for isomer 2..
As expected, the numbers in Fig. 7 have near bilateral symmetry with respect to a mirror plane drawn through C1 and C9. Values of \(\langle {\hat{S}}^2 \rangle\) vary from 0.00 to about 1.00 (in some cases we see \(\langle {\hat{S}}^2 \rangle\) slightly exceeding 1.00). This later value is typical of what is expected for a diradical. Interestingly, the sites with zero spin contamination seem almost, but not quite, to alternate with sites with significant spin contamination. They are also mostly on the front side, rather than the back side, relative to C1. The value of \(\langle {\hat{S}}^2 \rangle = 1\) has already been explained by the Lewis dot structure (LDS) shown in Fig. 6. The implications for further reactivity with a third \(^\bullet\)OH radical are evident. However, we should remember that the unpaired spins may be widely distributed around the molecule in a way that corresponds to multiple diverse resonance structures in the Lewis representation. The reader is invited to explore which different LDSs may be drawn for each isomer.

Figure 7: Spin contamination: Averaged over gas-phase and aqueous-phase values,the numbers show the degree of spin-contaminationat a C-O distance of 1.4 Å for the second OH group. Notethat the first OH group is always on carbon 1. It is important tonote that the underlying figure is the same as in Fig. 1because the carbon numbers are obscured by the red and blue numbersshowing the value of \(\langle{\hat{S}}^2 \rangle\)..
Inspection of the "BDE" at \(R\)(C-O) = 1.4 Å relative to a zero of energy at large \(R\)(C-O) for symmetry pairs showed that our results are plagued by more numerical problems than is the case for spin contamination. Often these are in cases where spin contamination is high and there might be more than one way to break symmetry. When the BDEs at \(R\)(C-O) = 1.4 Å for symmetry pairs are very different, we also usually find significant differences between the corresponding PECs and \(\langle {\hat{S}}^2 \rangle\) curves (shown in the SI). These numerical problems translate into an analysis problem that we have only been able to solve by “cleaning” (i.e., removing from further analysis) symmetry pairs whose BDEs at \(R\)(C-O) = 1.4 Å differ by more than 4 kcal/mol. Specifically, all of the symmetry pairs indicated by an asterisk in Table 1 have been removed from the data. The cleaned results for both \(\langle {\hat{S}}^2 \rangle\) and the BDE at \(R\)(C-O) = 1.4 Å are shown in Fig. 8. It is seen that \(\langle {\hat{S}}^2 \rangle\) is the same in gas phase and in aqueous phase. The gas-phase and aqueous-phase BDE are now linearly related with the BDE being 0.86 kcal/mol higher in the gas phase than in the aqueous phase. This is consistent with the idea that the products, \(^\bullet\)C\(_{60}\)(OH)\(_2\) and (especially) \(^\bullet\)OH are significantly stabilized relative to the reactant C\(_{60}\)(OH)\(_2\) in aqueous solution.

Figure 8: Aqueous-phase versus gas-phase correlation graph of\(\langle \hat{S}^2 \rangle\) and BDEs at \(R\)(C-O) = 1.4 Å..
Figure 9 provides a closer look at the correlation between the BDE and \(\langle {\hat{S}}^2 \rangle\) in both gas phase and in aqueous solution. For the isomers where \(\langle {\hat{S}}^2 \rangle=0\), there is a spread of BDEs governed by the chemical reactivity indices discussed above. For \(\langle {\hat{S}}^2 \rangle > 0\), there is a decrease in the BDE as the spin-contimination increases. The functional dependence of the BDE as a function of \(\langle {\hat{S}}^2 \rangle\) seems to be essentially the same in gas and in aqueous solution, even if the BDE in the gas phase is more highly bound by 0.86 kcal/mol.

Figure 9: BDEs as a function of \(\langle \hat{S}^2 \rangle\) at \(R\)(C-O) = 1.4 Åafter data cleaning..
By now it should be pretty clear that the PECs in gas-phase and in aqueous-phase for a given isomer are very similar. This is illustrated by Fig. 10 and further illustrated by figures for other isomers given in the SI.

Figure 10: Comparison of the short- and long-range shapes of the solution1,5-C\(_{60}\)(OH)\(_2\) (isomer 6) vertical dissociation PEC withthe gas-phase PEC: left, PECs matched at short and long range;right, corresponding solution minus gas phase PEC difference curves..
Lastly, in our exploratory work, we wanted to identify at least some of the responsible for the shape of the PECs. We expected the long-distance behavior to go as \(1/R^3\) which is typical of a dipole-dipole interaction. Instead we were surprised to see a region with a \(1/R^6\) van der Waals interaction-like behavior (Fig. 11). There is a simple argument that DFT should be unable to describe the van der Waals \(1/R^6\) part of the potential when the electron density of the separated fragments no longer overlap. Our calculations do not contradict this idea as the \(1/R^6\) is no longer present at very long distance. However, this is a reminder that a good functional can capture some of the \(1/R^6\) behavior at moderately long distance. Indeed one of the challenges when adding van der Waals corrections to DFT is not to overcount any \(1/R^6\) behavior that may be already present in the functional.

Figure 11: Identification of the van der Waals region of the 1,5-C\(_{60}\)(OH)\(_2\)(isomer 6) vertical dissociation PECs. The upper row shows log-logplots of the PEC and tries to “fit” the van der Waals region. The lowerrow uses a finite difference method to calculate the derivative and thendoes a least-squares fit to a cubic polynomial..
Buckminsterfullerene (C\(_{60}\)) is such an excellent radical scavanger that it has been called a “radical sponge” [1]. This makes it interesting as a potential antioxidant for biological applications provided that its solubility can be improved. One way to improve the solubility of C\(_{60}\) in water is to add hydroxyl groups which may be done in a variety of ways, including by successive addition of hydroxyl radicals (\(^\bullet\)OH). The relative stability of different fullerenols (\(^\bullet\))C\(_{60}\)(OH)\(_n\) (which is a radial for odd \(n\)) has been studied theoretically in previous gas-phase work [16]. Whereas the previous work treated additional fullerenols, the present work is limited to fullerendiols, C\(_{60}\)(OH)\(_2\). However the present work does extend the previous work in three ways.
Firstly, the SMD implicit solvent model has been used to include the dielectric nature of water, albeit without including the effects of hydrogen bonding. As C\(_{60}\)(OH)\(_2\) is nearly spherical, we expected that the SMD model would give essentially the same result as the Onsager model which is based upon a spherical cavity. This was shown to be the case. The effect of the implicit solvent model on the dipole moment of the molecule is large and yet calculated reactivity indices show exactly the same qualitative trends as in the gas phase, confirming that \(^\bullet\)OH is an electrophilic radical [16]. All of these calculations were done using a spin-restricted formalism for C\(_{60}\)(OH)\(_2\) and a spin-unrestricted formalism for \(^\bullet\)C\(_{60}\)(OH) and for \(^\bullet\)OH.
Secondly, we calculated potential energy curves (PECs) for the dissociation reaction C\(_{60}\)(OH)\(_2\) \(\rightarrow\) \(^\bullet\)C\(_{60}\)(OH) + \(^\bullet\)OH in both gas-phase and with the SMD implicit solvent model. This required us to use a spin-unrestricted formalism throughout the reaction path so that the spins could pair up in the same orbitals in C\(_{60}\)(OH)\(_2\) but separate into different orbitals in the products \(^\bullet\)C\(_{60}\)(OH) + \(^\bullet\)OH. That still left us with a tremendous problem of the dissociation pathway of 59 C\(_{60}\)(OH)\(_2\) isomers, each of which has multiple conformers. We decided to first carry out a crude vertical dissociation of each isomer which already gave us some very interesting information. In particular, the gas-phase and aqueous-phase PECs turn out to be quantitatively similar.
Thirdly, we noticed that, contrary to our (and we suppose also that of other workers in this field) expectations, it is not always true that \(\langle \hat{S}^2 \rangle = 0\) for C\(_{60}\)(OH)\(_2\). In fact, it is only true for about a quarter of the isomers. We were able to convince ourselves by examining spin densities and by drawing resonance structures that 1,8-C\(_{60}\)(OH)\(_2\) has diradicaloid character and we expect that there must also be di- or even pluriradicaloid character for many of the other fullerenediols. This actually creates a new problem that we did not attempt to address, namely that it is well established that there are several different ways to break symmetry in magnetic clusters and finding the lowest energy broken energy solution is nontrivial. What we did do was to recognize that all but 3 of the fullerenediol isomers come in symmetry pairs and that these almost always give quantitatively similar PECs and graphs of \(\langle \hat{S}^2 \rangle\) along the dissociation pathway. This allows us to identify and discard pathways where the symmetry pairs behave too differently (which was the case for only 4 of the isomers). None of the PECs show any convincing evidence of a barrier along their PEC.
While we find it reassuring that many of our conclusions drawn from gas-phase calculations are confirmed in our implicit solvent model, the most important conclusion of the present work is the realization that arbitrary fullerenols should and do have radicaloid character which should be taken into account by symmetry breaking (as we have done) or by yet more sophisticated linear combination of configuration state functions within some sort of wave function or hybrid wave function/DFT approach.
We gratefully acknowledge travel funding for AJE from the U.S.-Africa Initative [35]. This work was partially funded by NIH R01GM108583 and R35GM151951 to GAC. Computing time for this project was provided by the University of North Texas CASCaM high-performance clusters NSF Grant Numbers CHE1531468 and OAC-2117247, by the NSF Extreme Science and Engineering Discovery Environment, ACCESS project Number TG-CHE160044, and the University of Texas at Dallas’ Cyberinfrastructure and Research Services. AP would like to thank Pierre Girard for technical support in the context of the Grenoble Centre d’Experimentation du Calcul Intensif en Chimie (CECIC) computers used for the Orca calculations reported here.