An exceptional functionalization of doped fullerene observed via theoretical studies on the interactions of sulfur-doped fullerenes with halogens and halides

Abstract

This work reports the study on the interactions of sulfur-doped fullerenes with halogens and halides (except iodine and iodide) to obtain a deep understanding of such interactions for employing them in possible applications such as sensors and surface adsorption. The ωB97XD DFT code has been utilized in this study to obtain adsorption energies in the gas phase and solvent. The energy outcomes showed good results for adsorption in both gas and solvent (using the PCM model). However, the gas phase interactions are more thermodynamically favorable. The formation of halide complexes releases more energy than the formation of halogen complexes and the strongest interaction belongs to the interactions of disulfur-doped fullerenes with halides. Donor–acceptor transitions are mostly affected by sulfur doping, which made the C–S bond an auxiliary tool for the absorption process. Density of state (DOS) plots demonstrate the enhanced modification of conductivity properties upon sulfur-doping on fullerene structures. Electron densities and their Laplacians at bond critical points (BCPs) of interaction sites, NCI calculations and further visualizations clearly prove the existence of these interactions. It is worth noting that, in some cases, some observations like the partial functionalization of fullerene was seen and these observations were proven via QTAIM, NCI and energy data. In these cases, the stable and thermodynamically favorable cation of the halogenated doped fullerene (along with X₃⁻ as the counter ion) could be produced from the interaction of double-doped fullerene with two molecular halogens.

---

Introduction

Since fullerenes were first synthesized in 1985 by Smalley and Croto, new horizons have been expanded in the chemistry of carbon nanomaterials.¹ Compared to diamond and graphite, fullerenes are neutral, edgeless and have neither dangling bonds nor unpaired electrons. This gives fullerenes unique properties such as high-speed twirling around itself.² ARC discharge technique has been used as the main route for the synthesis of C₆₀ and C₇₀ fullerenes.^3,4 These molecules can be employed in superconductive materials,⁵ organic electric conductors,⁶ gas reservoirs, gas separators, fuel cells and lubricants.² Tegos and coworkers showed the photosynthetic and antimicrobial properties of functionalized fullerenes by selective photodynamic treatment to produce active oxygen.⁷

Although the mentioned applications of fullerenes are unique, the heterodoping of fullerenes gives them new potency in electronic instruments, nanocomposites and sensors.^8,9 Wang et al. reported a study on the adsorption process on C₃₅B–H₂–C₃₅B (B-doped fullerene) from two directions. They observed barrierless H₂ bond breaking from physisorption toward chemisorption via 3D-potential energy curves.¹⁰ In the case of sulfur-doping, Glenis et al. succeeded in the synthesis of C_60−2nS_n, C_60−3nS_n and C_60−4nS_n fullerenes from thiophene using graphite electrode ARC discharge.¹¹ Moreover, density functional studies were employed to calculate the energies and structures of thiafullerenes in several oxidation states. The results predicted that C₅₉S and C₅₉S²⁺ in the closed-cage form and C₅₉S⁴⁻ in the opened-cage structure are the major forms among the possible structures.¹² Sulfur-doping enhances the semiconductor behavior of CNTs, graphene and fullerene.^13–16 The adsorption of molecular iodine on sulfur-doped fullerenes (SFs) via noncovalent interactions has been reported using QTAIM and NBO analyses. The report showed considerable interactions between iodine and S-doped fullerene, which could be used in sensors.¹⁶ Therefore, doped fullerenes are suitable candidates to employ in the adsorption of various molecules and preparation of new sensors. The best tool to study these interactions is computational chemistry, which could precisely predict the existence and strength of noncovalent interactions as well as has the ability to predict chemical properties, geometric structures and compound specifications.¹⁷ Moreover, in these systems, using computational tools, three types of noncovalent interactions could be investigated, which include the halogen bond (between X–X and S with the X–X⋯S angle near to 180 degrees), π-halogen bond (between the halogen and π bond) and chalcogen bond (between C–S and halide or halogen, halogen with a 90 degree angle and the halogen is the electron donor).^18–20 Therefore, in continuation of our previous studies in computational chemistry,^21,22 especially on the study of noncovalent interactions,^14–16 we have decided to attempt a complete study on the interactions of SFs with halogens and halides using appropriate DFT calculations. Because, despite the fact that there are several studies on the adsorption or sensor potencies of simple and doped fullerenes (or other nanostructures),^23–26 by reviewing the literature, it was found that there are no reports related to the interaction of halogens or halide anions with doped fullerenes.

This study will be a guide to realize the interactions between SFs and halogens or halides (except iodine and iodide) and propose possible applications that may be used in the future for several purposes such as drug delivery, preparation of sensors, and storage. Moreover, QTAIM, NBO and NCI calculations were employed in addition to DFT calculations to provide more evidence about the nature and strength of these interactions. Moreover, we provide some evidence related to the functionalization (with halogens) of SFs in some special cases. Details of the computations and results obtained in this work are presented in the following sections.

Methods

To start the simulation of noncovalent interactions between SFs and halogens or halides, we optimized C₆₀ fullerene, and after geometry optimization, one and two sulfur atoms were substituted permanently to obtain C₅₉S (SF) and C₅₈S₂ fullerenes (S2F1 and S2F2), respectively. The next step was to place molecular fluorine, chlorine and bromine (with the code name X₂ (X = F, Cl and Br)) and halide anions (fluoride, chloride and bromide as X (X = F⁻, Cl⁻ and Br⁻)) on the surfaces of the doped fullerenes and optimize their structures. It should be noted that in monosulfur-doping, the model name was set to SF and in the disulfur-doped fullerenes, sulfur was located approximately 90° and 180° from each other toward the center of the fullerene sphere, which were named as S2F1 and S2F2, respectively. Density functional theory (DFT) method was employed in this work because it is vastly used to study noncovalent interactions due to its lower complexity against an increasing system size, which is an advantage of DFT method. Moreover, DFT calculations have been shown to have high abilities in the calculation of molecular properties^27,28 of chemical structures, which are comparable with the most computationally expensive MP2 methods,^29,30 and DFT calculations can also be used for the prediction of thermodynamic and kinetic properties of chemical phenomena.^31,32 Among the various DFT methods, ωB97XD, designed by Chai and Head-Gordon,^33,34 is the best DFT code to interpret intermolecular interactions because many DFT methods cannot describe dispersive energy. This method is a development of the former ωB97X method, obtained by applying some empirical corrections for binary atomic dispersion.^35,36 After optimizations, the interaction energies for all the complexes were calculated using eqn (1).

ΔE_ads = E_complex − (E_absorbent + E_adsorbate) (1)

Furthermore, to investigate solvent effects, PCM model was used to obtain solvation energies for every single optimized complex in benzene, chloroform and cyclohexane (three solvents with different polarities).³⁷ Eqn (1) was employed to determine the interaction energies in the solvents. All the geometry optimizations and calculation of the solvent effects were performed using the Gaussian 09 Revision A.01 suite³⁸ with the ωB97X-D/6-31G level of theory. We used this very small basis set to decrease the calculation costs (because of the use of large systems), and more importantly we have previously shown that there is no meaningful difference in the energy results when the level of the basis set is increased in DFT calculations.³⁹

With regards to the valuable work of Weinhold et al. for the introduction of NBO concepts, this calculation was also performed on each optimized structure to yield accurate atomic charges and donor–acceptor transactions of SFs with halogens and halides via NBO 3.1, integrated in the Gaussian 09 program.⁴⁰ To monitor the changes in conductive levels of SFs, density of states (DOS) plots were depicted by utilizing the GaussSum software.⁴¹ To calculate reactivity parameters, Koopman's theorem⁴² was employed to obtain the chemical potential (μ), chemical hardness (η), global softness (S) and electrophilicity index (ω) for all the complexes from the following equations.

μ = (E_LUMO + E_HOMO)/2 (2)

η = (E_LUMO − E_HOMO)/2 (3)

S = 1/η (4)

ω = μ²/2η (5)

Bader's quantum theory of atoms in molecules (QTAIM)⁴³ was employed to calculate the electron densities and Laplacians of electron densities at bond critical points (BCPs) in noncovalent interaction sites using the AIMAll program.⁴⁴ Finally, based on the interesting studies of Johnson, Contreras-Garcia and colleagues on noncovalent interaction (NCI) indexes, it seems strongly indispensable to perform this calculation in the current research as well. The NCI indexes are related to the fundamental dimensionless parameter in DFT, which is called the reduced density gradient (RDG), as described by the S parameter in eqn (6).

Since the sign of λ₂ is proportional to the essence of the interactions, the plots consisting of the sign of λ₂ × ρ versus RDG would indicate a noncovalent interaction near zero area (horizontal axis).⁴⁵ To fulfil this, the wavefunctions obtained from the optimized structure for each complex was examined using the NCIPLOT program⁴⁶ and depicted using the VMD 1.9.1 software,⁴⁷ in which the green color in the isosurfaces means a weak and noncovalent interaction between the SFs and halogen or halides. The plots of sign λ₂ × ρ versus RDG were drawn using gnuplot 4.6.5.⁴⁸

Result and discussion

Optimized structures

The first step of this study is the optimization of three doped fullerenes (SF, S2F1 and S2F2) to consider their interactions with halogens and halides. The optimized structures of these three initial structures (or adsorbents, generally named SFs) are shown in Fig. 1. Then, they were interacted with six molecules (three halides and three molecular halogens) as adsorbates, i.e. F⁻, Cl⁻, Br⁻, F₂, Cl₂ and Br₂. For the complexes with molecular halogens, it is possible to consider two different situations of halogens versus SFs, one is direct (the angle of SF–X–X is nearly 180 degrees, which could be a halogen bond) and another is perpendicular (the angle of SF–X–X is nearly 90 degrees, which could be a π-halogen or chalcogen bond). It is noticeable that in the complexes with halide anions, mostly the chalcogen bond exists.

Fig. 1 The optimized structures of mono and disulfur-doped fullerene as employed adsorbents in the calculations.

Therefore, for the complex between SF and each molecular halogen, two different configurations could be considered i.e., the direct position or SF–D–X₂ and perpendicular position or SF–P–X₂. Although, for the complexes of S2F1 (or S2F2) with each molecular halogen, three different configurations could be considered i.e., the direct position of both halogens or S2F1–2D–X₂, perpendicular position of both halogens or S2F1–2P–X₂ and direct position of one halogen with perpendicular position of another halogen or S2F1–DP–X₂. Therefore, a total of 33 complexes between the SFs and halogens or halides were considered and their structures were optimized. The optimized structures of these complexes are depicted in Fig. 2–4, for the SF, S2F1 and S2F2 complexes, respectively. In these figures, the minimum distance between the absorbent and adsorbate has been mentioned for each complex.

Fig. 2 The optimized structures for the noncovalent complexes of monosulfur-doped fullerene (SF) with halogens and halides.

Fig. 3 The optimized structures for the noncovalent complexes of S2F1 (first disulfur-doped fullerene) with halogens and halides.

Fig. 4 The optimized structures for the noncovalent complexes of S2F2 (another disulfur-doped fullerene) with halogens and halides.

The average C–S bond lengths in the SFs (obtained from the calculations) are 1.864 Å and 1.863 Å, for single and double doping, respectively. In the complexes with halogens, the direct positions with the halogen bond have the shorter SF⋯X distance versus the perpendicular positions with chalcogen or π-halogen bonds.

The expected structural deviation from the normal spherical form of fullerene was seen in the doped sulfur sites. Noticeably, it was found that despite the initial situation of each model, some complexes changed from the direct mode to the perpendicular situation or the reverse after the optimization process. Examples for this change are SF–D–Br₂, S2F2–2D–Cl₂, S2F2–2D–Br₂ and all S2F1–DP and S2F2–DP models. Two later cases showed unusual geometrical parameters such as small interatomic distances between one halogen atom and carbon atom that is next to the dopant sulfur. These cases could be considered as partial functionalization and will be discussed from the energy and QTAIM point of view in the following sections. Apparently, all the minimum distances between absorbent and adsorbates (except the functionalization cases) were between 1.8–3.7 Å. Therefore, these interactions could be classified as noncovalent interaction, according to the definitions given by Cerny and Hobza considering the distance of the noncovalent interaction.⁴⁹ As a general comparison, it should be mentioned that typically these distances are the lowest in magnitude for fluorine-bearing complexes. In the DP models, the lowest distances were observed in the S2F1–DP models, in which one halogen atom was located in the covalent distance of the carbon atom that was next to the doped sulfur atom. Furthermore, other data help us to explain such occurrence.

Interactions energies in gas and solvent

To determine the possibility and strength of the interactions, the adsorption energies of halogens and halides on the surfaces of SFs for various possible configurations of each model in the gas phase and three solvents (benzene, chloroform and cyclohexane) were calculated according to eqn (1) and the results are provided in Table 1.

Table 1 Interaction energies of all the complexes in gas and solvents (benzene, chloroform and cyclohexane)

ΔEads	Gas	Benzene	Chloroform	Cyclohexane
SF–F	−77.31	−46.40	−34.32	−49.39
SF–Cl	−28.68	−10.42	−4.22	−12.06
SF–Br	−32.81	−14.97	−8.74	−16.59
SF–D–F2	−1.14	−1.02	−0.95	−1.04
SF–D–Cl2	−1.15	−1.04	−0.95	−1.06
SF–D–Br2	−7.48	−7.46	−7.37	−7.47
SF–P–F2	−2.41	−2.31	−2.23	−2.32
SF–P–Cl2	−2.19	−2.07	−1.99	−2.09
SF–P–Br2	−7.43	−7.21	−7.08	−7.23
S2F1–2F	−93.99	−64.01	−52.85	−66.84
S2F1–2Cl	−8.54	0.51	2.14	−0.10
S2F1–2Br	−14.06	−7.48	−6.93	−7.85
S2F1–2D–F2	−0.30	−0.02	0.07	−0.04
S2F1–2D–Cl2	−21.76	−30.04	−35.21	−28.99
S2F1–2D–Br2	−26.97	−33.16	−37.22	−32.36
S2F1–2P–F2	−5.01	−4.82	−4.69	−4.84
S2F1–2P–Cl2	−6.42	−6.21	−6.09	−6.24
S2F1–2P–Br2	−20.77	−20.18	−19.77	−20.26
S2F1–DP–F2	−94.13	−97.38	−99.47	−96.96
S2F1–DP–Cl2	−52.87	−60.79	−65.27	−59.85
S2F1–DP–Br2	−62.82	−68.41	−71.69	−67.73
S2F2–2F	−104.56	−73.02	−60.75	−76.07
S2F2–2Cl	−14.39	−4.54	−2.24	−5.28
S2F2–2Br	−20.91	−12.41	−10.50	−13.04
S2F2–2D–F2	−2.79	−2.72	−2.67	−2.74
S2F2–2D–Cl2	−5.62	−5.60	−5.54	−5.61
S2F2–2D–Br2	−15.21	−14.73	−14.56	−14.76
S2F2–2P–F2	−9.92	−14.73	−17.85	−14.11
S2F2–2P–Cl2	−5.69	−5.31	−5.12	−5.34
S2F2–2P–Br2	−16.22	−15.88	−15.63	−15.93
S2F2–DP–F2	−22.62	−33.61	−40.56	−32.22
S2F2–DP–Cl2	−12.30	−16.42	−19.21	−15.88
S2F2–DP–Br2	−22.53	−26.21	−28.72	−25.72

---

The calculated results showed thermodynamically desirable energies for all noncovalent interactions in both the gas and solvent phase. Generally, the adsorption of halides is more favorable than molecular halogens by noticeable values. For the interactions with a halide, more noncovalent interaction energy was released with the order of fluoride > bromide > chloride and the maximum interaction energy was calculated for the S2F2–2F model with 104.56 kcal mol⁻¹. However, the double-doped complexes have adsorption energies slightly less than twice the single-doped complexes. This means that the relation between adsorption energy and the number of dopants is not completely linear. In the solvent, the highest stability belonged to the fluoride complexes, and in all the halide complexes, cyclohexane is the solvent with the highest released interaction energies. This could be assigned to the lower polarity of the complexes when they interact together versus when the complex parts remain individual in the solvent. Therefore, the interaction of each halide typically was more powerful than that of the molecular halogen analogues. It could be concluded that the chalcogen bond (exists in halide complexes) has greater strength than the halogen bond (exists in direct halogen complexes) and π-halogen bond.

For the interactions of halogens with SFs, we studied the SF–D and SF–P models. As long as the bromine molecule remained perpendicular to the surface of the fullerene, it could be deduced that interactions in which the halogen molecule stands perpendicularly above the fullerene have higher stability than the directly situated ones and is correct generally for the S2F1–2P and S2F1–2D models as well. In other words, bromine prefers the chalcogen or π-halogen bond rather than the halogen bond. For chlorine and bromine, noncovalent interactions with S2F1 are more favorable than S2F2, but for fluorine case this is reverse. Moreover, in complexes with one molecular halogen, perpendicular positions have more negative interaction energies, which show preference for π-halogen or chalcogen bond rather than the halogen bond in these systems.

More importantly, an interesting observation was observed for the S2F1–DP model in which a large number of stabilization energies were calculated, and it was not acceptable to sort them out according to noncovalent interaction classification. In these complexes, the energy values up to 3 times larger than the other analogues models, higher interaction energies in the solvent (opposite to the observed orders in other complexes) and small distances between the halogen atom and carbon atom adjacent to the doped sulfur suggest a new phenomenon that is beyond the simple noncovalent interactions. In these cases, it appears that the SFs were functionalized by halogens and we obtained [SF–X]⁺X₃⁻. This type of functionalization could be interesting and useful in future research with regards to the chemistry of doped fullerenes. However, more quantum mechanics-based analyses are needed to confirm this observation and we wish to provide them in the next sections.

NBO analyses

NBO calculations were executed to obtain some useful information about the interactions. Two main categories of data were obtained from the NBO calculations, which are atomic charges and donor–acceptor transaction (E₂ interaction energies). Table 2 charts the selected important atomic charges involved in the noncovalent interactions of SF complexes (only for fullerene doped with one sulfur atom) along with SF alone; moreover, the full results of NBO atomic charges could be observed in the ESI (Table S1†). In Table 2, the average of atomic charges of three carbon atoms connected to the doped sulfur atom was reported as the C (Av) charge, the charge of the doped sulfur atom was reported as the S charge and the charges of the halogen atom(s) (one atom for halide complexes and two atoms for halogen complexes) were reported as X(1) (the closer halogen atom) and X(2).

Table 2 Selected NBO atomic charges for complexes of SF with halogens and halides

Complexes	S	C^a (Av)	X(1)	X(2)
a This value is the average of the atomic charges of three carbon atoms connected to the doped sulfur.
SF (alone)	0.859	−0.198	—	—
SF–F	0.977	−0.178	−0.731	—
SF–Cl	0.934	−0.171	−0.866	—
SF–Br	0.918	−0.173	−0.820	—
SF–D–F2	0.871	−0.197	−0.026	−0.010
SF–D–Cl2	0.853	−0.194	−0.028	−0.004
SF–D–Br2	0.856	−0.196	−0.007	0.018
SF–P–F2	0.860	−0.197	−0.003	0.006
SF–P–Cl2	0.857	−0.196	−0.005	0.007
SF–P–Br2	0.852	−0.197	0.000	0.013

---

According to the atomic charge data, the negative charges (for the halogen atom) in all noncovalent interactions are in the order: chloride > bromide > fluoride. This means that fluoride transfers more negative charge to SFs versus the other halides. These values are in accordance with the calculated adsorption energies, which showed that fluoride has the best interaction (and the maximum charge transfer) and chloride has the worst interaction (and the minimum charge transfer) in our models. Moreover, the alteration of the atomic charge of the SFs (versus SF alone) in the halide complexes is more than that in the halogen complexes. For example, the S charge in SF, SF–F, SF–D–F₂ and SF–P–F₂ are 0.859, 0.977, 0.871 and 0.860 au, respectively. The same conditions could be observed for the C (Av) charges.

In the complexes with halogens (molecular), for the halogens connected directly to the fullerene, much more atomic charge was observed in contrast with those remaining perpendicular. This could be seen in the SF–D–X₂ and SF–P–X₂ (X = F, Cl) models, while in the bromine complexes, they are perpendicular in both modes and they had more positive charges compared with the fluorine and chlorine complexes. It is noticeable that for S2F1–DP–X₂ (X = F, Cl, Br) and S2F2–DP–X₂ (to a lesser extent), an inharmonic rhythm was observed, especially for the atomic charge of the carbon atoms that were next to the doped sulfur atom. Therefore, more negative charges were observed for the halogen atoms in contrast with the other 2D or 2P analogue models.

Table 3 presents the strongest second order perturbation energies for the donor–acceptor transactions of each model. In the halide complexes, in agreement with previous results, fluoride has the highest E₂ interaction energies, while in the halogen complexes, there is no meaningful relation between the type of halogen and E₂ energies. In most of the interactions, the lone-pair of halogen acts as the electron donor and the C–S antibonding sigma bond (σ* C–S) acts as the electron acceptor, which could be a reference to the role of S-doping in driving the transactions in such a manner. Moreover, except in some cases, the π-bond does not participate as a donor or acceptor in these interactions. These observations show that in the perpendicular positions, the chalcogen bond plays a major role in these interactions. However, in models that consist of the direct configuration of halogens versus SFs (except for the bromine complexes that stay perpendicular in both), the transaction could be observed from the lone-pair of doped-sulfur to the X–X antibonding electrons. Noticeably, some complexes, such as S2F1–2D–Br₂ and S2F2–DP–X₂ (X = F, Cl, Br), show large E₂ energies in their electron transfer from the lone-pair of the carbon atom (next to dopant sulfur) to the σ* halogens. In few cases (S2F1–2D–Cl₂ and S2F2–2P–F₂), one halogen atom is so close to the carbon that it makes them a unified acceptor group and also very strong transitions were observed for them. Finally, for the S2F1–DP–X2 (X = F, Cl, Br) models, the LP_X → σ* C–S transitions have higher energy values in contrast with similar transitions in the other models.

Table 3 The strongest second order perturbation energies (E₂) of donor–acceptor transactions for selected complexes

Complex	Donor	Acceptor	E2^a
a All values are in kcal mol^−1.
SF–F	LP F	σ* C–S	37.09
SF–Cl	LP Cl	σ* C–S	9.01
SF–Br	LP Br	σ* C–S	8.82
SF–D–F2	LP S	σ* F–F	2.97
SF–D–Cl2	LP S	σ* Cl–Cl	3.30
SF–D–Br2	LP Br	σ* C–S	2.45
SF–P–F2	LP F	σ* C–S	0.65
SF–P–Cl2	LP Cl	σ* C–S	0.87
SF–P–Br2	LP Br	σ* C–S	2.89
S2F1–2F	LP F	σ* C–S	35.71
S2F1–2Cl	LP Cl	π* C–C	10.78
S2F1–2Br	LP Br	σ* C–S	22.28
S2F1–2D–F2	LP F	σ* C–S	0.48
S2F1–2D–Cl2	LP Cl	σ* C–Cl	117.57
S2F1–2D–Br2	LP C	σ* Br–Br	239.28
S2F1–2P–F2	LP F	σ* C–S	0.95
S2F1–2P–Cl2	LP Cl	σ* C–S	0.40
S2F1–2P–Br2	LP C	σ* Br–Br	4.53
S2F1–DP–F2	LP F	σ* C–S	45.34
S2F1–DP–Cl2	LP Cl	σ* C–S	16.41
S2F1–DP–Br2	LP Br	σ* C–S	16.17
S2F2–2F	LP F	σ* C–S	64.40
S2F2–2Cl	LP Cl	σ* C–S	16.20
S2F2–2Br	LP Br	σ* C–S	19.21
S2F2–2D–F2	LP S	σ* F–F	3.47
S2F2–2D–Cl2	LP Cl	π* C–C	0.25
S2F2–2D–Br2	LP Br	σ* C–S	2.48
S2F2–2P–F2	LP F	σ* C–F	123.56
S2F2–2P–Cl2	LP Cl	σ* C–S	0.93
S2F2–2P–Br2	LP Br	σ* C–S	3.35
S2F2–DP–F2	LP C	σ* F–F	204.29
S2F2–DP–Cl2	LP C	σ* Cl–Cl	172.70
S2F2–DP–Br2	LP C	σ* Br–Br	165.70

---

Population analyses and DOS plots

To realize the sensor behavior of SFs, by probing their HOMO and LUMO energies and energy gaps, it appears essential to obtain some reactivity indexes from eqn (2) to (5) to describe the change in chemical behavior upon S-doping. Table 4 provides the HOMO and LUMO energy values, their energy gaps, chemical potential, chemical hardness, chemical softness and electrophilicity indexes of fullerene (C₆₀) and employs SFs to compare the behavior of fullerene after sulfur-doping. In addition, to obtain more insight about the electronic properties of these structures, Fig. 5 shows the DOS plots for each of them (C₆₀ and SFs).

Table 4 Energies of HOMO, LUMO, energy gap (E_g), chemical potential (μ), chemical hardness (η), global softness (S) and electrophilicity index (ω) for fullerene and SFs. All the energy values (E) are in the eV scale

Fullerene	EHOMO	ELUMO	Eg	μ	η	S	ω
F (C60)	−0.2984	−0.0713	0.2270	−0.1849	0.1135	8.80	0.1505
SF	−0.2571	−0.0714	0.1856	−0.1643	0.0928	10.77	0.1454
S2F1	−0.2519	−0.0722	0.1797	−0.1620	0.0898	11.13	0.1461
S2F2	−0.2554	−0.0716	0.1838	−0.1635	0.0919	10.88	0.1455

---

Fig. 5 DOS plots of fullerene and SFs (HOMO and LUMO levels are indicated by the narrow bar). Energy gap dimensions are in eV.

As determined, C₆₀ fullerenes have a high E_g value versus SFs, and upon doping their E_g value is reduced; this decrease is larger when the number of adjacent dopant atoms increases. Moreover, the HOMO energy level shifts up (to more positive values) upon doping, which means an increase in nucleophilicity or electron donation property and the LUMO energy is reduced a little (to more negative values); this means that there is a decrease in the electrophilicity index (this value is confirmed by the electrophilicity index in the last column in Table 4).

These values show the reduced reactivity (and increase in stability) of fullerene after doping, which is confirmed by the chemical potential (μ) values listed in the 5^th column in Table 4. The observed E_g values for SF and S2F2 are similar to each other but for the S2F1 fullerene, a lower E_g gap was revealed. The obtained results show the modification of the electronic properties due to the sulfur doping in fullerenes from about 18–20%. Surely, the chemical hardness decreases and global softness increases by the sulfur doping of C₆₀ fullerene. All these values confirm the distinct change in the behavior of SFs versus simple fullerenes and these can be used in the desired manner to obtain better results.

QTAIM analyses and NCI indexes

As mentioned before, the electron density (ρ) and its Laplacian (∇²ρ) at the BCPs of noncovalent interaction sites could ensure the existence of this type of interactions and also determine their characteristics. After performing the QTAIM calculation, ρ and ∇²ρ values from the BCPs, which were desired on interaction area, were extracted. Table 5 presents the important data of electron densities and their Laplacians at the BCPs in the noncovalent interaction sites between SFs and halogens or halides. In some cases, more than one critical point was observed for each complex. In addition, a unique picture of the noncovalent interaction for each case with clarifying green and red dots on the pictures, which represent the BCPs and RCPs (Ring Critical Point), are shown in the ESI (Fig. S1–S3†).

Table 5 Electron densities (ρ) and their Laplacian (∇²ρ) of important BCPs at noncovalent interaction sites for all complexes obtained from the QTAIM calculations

∇^2ρ(e/ao^5)	ρ(e/ao^3)	Type	Complex
S2F1–2F	F–S	5.5 × 10^−2	1.8 × 10^−1
	F–S	5.4 × 10^−2	1.8 × 10^−1
S2F1–2Cl	Cl–C	1.9 × 10^−2	5.4 × 10^−2
	Cl–S	1.6 × 10^−2	5.3 × 10^−2
S2F1–2Br	Br–S	2.5 × 10^−2	6.7 × 10^−2
	Br–S	2.2 × 10^−2	6.0 × 10^−2
S2F2–2D–F2	F–S	1.9 × 10^−2	7.7 × 10^−2
	F–S	1.2 × 10^−2	7.2 × 10^−2
S2F2–2D–Cl2	Cl–C	4.6 × 10^−3	1.5 × 10^−2
	Cl–C	4.1 × 10^−3	1.4 × 10^−2
	Cl–C	4.6 × 10^−3	1.5 × 10^−2
	Cl–C	4.3 × 10^−3	1.4 × 10^−2
S2F2–2D–Br2	Br–S	9.6 × 10^−3	3.3 × 10^−2
	Br–C	8.1 × 10^−3	2.5 × 10^−2
	Br–C	6.2 × 10^−3	1.9 × 10^−2
	Br–S	9.5 × 10^−3	3.2 × 10^−2
	Br–C	8.2 × 10^−3	2.5 × 10^−2
	Br–C	6.0 × 10^−3	1.9 × 10^−2
S2F2–2P–F2	F–C	9.3 × 10^−2	2.5 × 10^−1
	F–S	6.1 × 10^−3	3.0 × 10^−2
	F–C	5.3 × 10^−3	2.4 × 10^−2
S2F2–2P–Cl2	Cl–C	5.0 × 10^−3	1.6 × 10^−2
	Cl–S	4.7 × 10^−3	1.8 × 10^−2
	Cl–C	4.5 × 10^−3	1.5 × 10^−2
	Cl–C	4.5 × 10^−3	1.5 × 10^−2
S2F2–2P–Br2	Br–S	1.1 × 10^−2	3.9 × 10^−2
	Br–S	7.7 × 10^−3	2.8 × 10^−2
	Br–C	7.5 × 10^−3	2.8 × 10^−2
	Br–S	9.3 × 10^−3	3.2 × 10^−2
	Br–C	8.0 × 10^−3	2.5 × 10^−2
S2F2–DP–F2	F–C	1.3 × 10^−1	2.9 × 10^−1
	F–C	7.0 × 10^−2	2.2 × 10^−1
S2F2–DP–Cl2	Cl–C	5.7 × 10^−2	1.1 × 10^−1
	Cl–S	4.9 × 10^−3	1.9 × 10^−2
	Cl–C	4.3 × 10^−3	1.4 × 10^−2
	Cl–C	3.3 × 10^−3	1.1 × 10^−2
S2F2–DP–Br2	Br–C	6.3 × 10^−2	9.6 × 10^−2
	Br–S	1.1 × 10^−2	4.0 × 10^−2
	Br–C	7.5 × 10^−3	2.2 × 10^−2
	Br–S	7.4 × 10^−3	2.7 × 10^−2
S2F2–2F	F–S	7.7 × 10^−2	2.3 × 10^−1
	F–S	5.5 × 10^−2	1.8 × 10^−1
S2F2–2Cl	Cl–S	2.2 × 10^−2	6.6 × 10^−2
	Cl–S	1.6 × 10^−2	5.2 × 10^−2
	Cl–C	8.9 × 10^−3	2.7 × 10^−2
S2F2–2Br	Br–S	2.3 × 10^−2	6.3 × 10^−2
	Br–C	2.2 × 10^−2	6.0 × 10^−2
SF–F	F–S	5.7 × 10^−2	1.9 × 10^−1
	F–C	2.3 × 10^−2	9.8 × 10^−2
	F–C	2.3 × 10^−2	9.8 × 10^−2
SF–Cl	Cl–S	1.9 × 10^−2	5.8 × 10^−2
	Cl–C	1.2 × 10^−2	4.0 × 10^−2
SF–Br	Br–S	1.9 × 10^−2	5.5 × 10^−2
	Br–C	1.5 × 10^−2	4.6 × 10^−2
SF–D–F2	F–S	1.8 × 10^−2	4.7 × 10^−2
SF–D–Cl2	Cl–S	1.2 × 10^−2	1.5 × 10^−2
SF–D–Br2	Br–S	9.6 × 10^−3	3.3 × 10^−2
	Br–C	8.0 × 10^−3	2.5 × 10^−2
	Br–C	6.3 × 10^−3	2.0 × 10^−2
SF–P–F2	F–C	5.5 × 10^−3	2.4 × 10^−2
	F–S	4.8 × 10^−3	2.7 × 10^−2
	F–S	4.1 × 10^−3	2.3 × 10^−2
SF–P–Cl2	Cl–S	4.9 × 10^−3	1.9 × 10^−2
	Cl–C	4.5 × 10^−3	1.4 × 10^−2
SF–P–Br2	Br–S	1.1 × 10^−2	3.8 × 10^−2
	H–C	8.1 × 10^−3	2.7 × 10^−2
S2F1–2D–F2	F–S	4.7 × 10^−3	2.3 × 10^−2
	F–S	4.2 × 10^−3	2.2 × 10^−2
S2F1–2D–Cl2	Cl–C	9.3 × 10^−2	1.1 × 10^−1
	Cl–C	1.4 × 10^−2	5.0 × 10^−2
S2F1–2D–Br2	Br–C	8.7 × 10^−2	8.8 × 10^−2
	Br–S	1.3 × 10^−2	4.8 × 10^−2
S2F1–2P–F2	F–S	5.9 × 10^−3	3.0 × 10^−2
	F–C	4.1 × 10^−3	1.9 × 10^−2
	F–C	4.0 × 10^−3	1.9 × 10^−2
	F–C	4.7 × 10^−3	2.2 × 10^−2
	F–C	4.7 × 10^−3	2.2 × 10^−2
	F–S	4.4 × 10^−3	2.4 × 10^−2
S2F1–2P–Cl2	Cl–C	5.2 × 10^−3	1.6 × 10^−2
	Cl–C	4.4 × 10^−3	1.4 × 10^−2
	Cl–C	4.8 × 10^−3	1.5 × 10^−2
	Cl–C	4.6 × 10^−3	1.5 × 10^−2
S2F1–2P–Br2	Br–C	1.2 × 10^−2	3.7 × 10^−2
	Br–S	1.2 × 10^−2	4.2 × 10^−2
	Br–C	1.1 × 10^−2	3.3 × 10^−2
	Br–C	1.1 × 10^−2	3.5 × 10^−2
	Br–S	9.8 × 10^−3	3.5 × 10^−2
	Br–S	7.1 × 10^−3	2.5 × 10^−2
	Br–C	6.7 × 10^−3	1.9 × 10^−2
S2F1–DP–F2	F–C	2.3 × 10^−1	−2.8 × 10^−1
	F–S	7.1 × 10^−2	2.3 × 10^−1
S2F1–DP–Cl2	Cl–C	1.5 × 10^−1	−5.8 × 10^−2
	Cl–S	2.4 × 10^−2	7.3 × 10^−2
S2F1–DP–Br2	Br–C	1.3 × 10^−1	−4.5 × 10^−2
	Br–S	2.4 × 10^−2	6.7 × 10^−2
	Br–S	1.2 × 10^−2	4.3 × 10^−2

---

It was found that electron densities in the critical points of the halide complexes are much higher than those in the halogen complexes. With regards to this, for the fluoride complexes, F–S and F–C interactions have the highest ρ values among the observed BCPs in the SF and S2F2 complexes. In the chloride cases, the interactions labeled with Cl–S and Cl–C have the maximum ρ values in the S2F1 and S2F2 complexes, respectively. However, Br–S and Br–C interactions have the highest ρ values in their BCPs in the S2F2 and SF complexes. The high magnitude and positive sign of the electron density's Laplacian in the BCPs indicates the electrostatic (noncovalent) essence of these interactions. For the S2F1–DP–X₂ (X = F, Cl, Br) complexes, which are described separately, the interaction of molecular fluorine with S2F2 has the highest ρ values in both F–S and F–C interactions.

Moreover, in SF and S2F1 complexes, Cl–S and Cl–C interactions, respectively, have the largest ρ values in their BCPs for the chlorine complexes. Finally, for bromine complexes, both Br–S and Br–C have the maximum magnitude of ρ in the BCPs by interaction with S2F1. All the above mentioned complexes have positive Laplacian of ρ at the BCPs, which illustrate that electrostatic interactions occur through this procedure.

However, something special was observed in S2F1–DP–X₂ (X = F, Cl, Br) complexes, in which some ρ values and the sign of their Laplacian suggested that it is not a formal noncovalent interaction. The ρ values for the X–C (X = F, Cl, Br) interactions were 0.2325, 0.1492, and 0.1338 e/a_o³ for fluorine, chlorine and bromine, respectively. Apparently, these values are not in the range of other complexes that bear electrostatic noncovalent interactions. This was confirmed when we observed a negative sign for the Laplacian of ρ value in relative BCPs which is −0.2841, −0.0576 and −0.0451 e/a_o⁵ for fluorine, chlorine and bromine, respectively. Table 6 presents an extended version of all BCP near interaction sites for all S2F1–DP complexes to give additional insight for understanding this phenomenon. As can be seen from the abovementioned table, the carbon next to the doped sulfur atom appears to have a covalent bond with one halogen atom due to the negative Laplacian of ρ at the given BCP. According to these data, a comparison between the ρ values of C–S BCPs indicates the reduced covalent strength of one C–S bond. With respect to this fact, since two halogen atoms interact with SF either covalently or noncovalently, the two remaining halogen atoms have some type of secondary interaction that is notable in magnitude (for F–F, Cl–Cl and likely Br–Br). Therefore, the low distance between halogens and SF will lead to partial functionalization (halogenation) of S2F1 and we obtain the stable structure as [S2F1–X]⁺X₃⁻, which is produced from S2F1 and two molecular halogens via an exothermic process. A more precise look at the results of the diminished C–S BCP's ρ values for fluorine, chlorine and bromine reveals the fact that fluorine would reduce the strength of the C–S bond more than chlorine and bromine.

Table 6 Electron densities (ρ) and their Laplacian (∇²ρ) of important BCPs for partially functionalized S2F1–DP–X₂ (X = F, Cl, Br) complexes

S2F1–DP–F2			S2F1–DP–Cl2			S2F1–DP–Br2
Type	ρ(e/ao^3)	∇^2ρ(e/ao^5)	Type	ρ(e/ao^3)	∇^2ρ(e/ao^5)	Type	ρ(e/ao^3)	∇^2ρ(e/ao^5)
F–F	0.1792	0.7514	Cl–Cl	0.0415	0.1125	Br–Br	0.0411	0.0923
F–F	0.0492	0.2661	Cl–Cl	0.0663	0.1195	Br–Br	0.0496	0.0885
F–F	0.0232	0.1065	Cl–Cl	0.0060	0.0212	Br–Br	0.0172	0.0559
C–F	0.2355	−0.2841	C–Cl	0.1492	−0.0576	S–Br	0.0236	0.0666
S–F	0.0713	0.2306	S–Cl	0.0243	0.0730	S–Br	0.0121	0.0425
S–C	0.1498	−0.1520	S–C	0.1472	−0.1332	C–Br	0.1338	−0.0451
S–C	0.1582	−0.1765	S–C	0.1587	−0.1746	S–C	0.1482	−0.1357
S–C	0.1421	−0.1067	S–C	0.1468	−0.1362	S–C	0.1593	−0.1754
S–C	0.1547	−0.1545	S–C	0.1501	−0.1387	S–C	0.1452	−0.1313
S–C	0.1095	−0.0318	S–C	0.1194	−0.0516	S–C	0.1490	−0.1358
S–C	0.1615	−0.1669	S–C	0.1566	−0.1473	S–C	0.1227	−0.0579
						S–C	0.1547	−0.1420

---

In the last step of this study, NCI index calculations were employed to produce noncovalent isosurfaces between absorbent and adsorbate species, which gave additional evidence for the noncovalent interactions in our complexes. These type of calculations have been recently developed as another proof for noncovalent interactions, which sometimes produce different results from the QTAIM calculations.⁵⁰ The diagrams of the isosurfaces (from several views for each model) for noncovalent interactions in the S2F1 complexes are shown in Fig. 6, and for the SF and S2F2 complexes, they have been shown in the ESI (Fig. S4 and S5†). In addition, Fig. 7 presents plots of sign λ₂ × ρ versus reduced density gradient (RDG) for all the complexes.

Fig. 6 Noncovalent interaction isosurfaces obtained from RDG and electron density frames for the interactions of the S2F1 model with halogens and halides.

Fig. 7 Plots of sign λ₂ × ρ versus RDG of the noncovalent interaction for all complexes.

In the abovementioned figure, the green-colored isosurfaces clearly show the noncovalent interaction by subtraction of the RDG frames from the electron density frames. According to the NCIPLOT reference,⁴⁵ green isosurfaces represent a weak van der Waals interaction between both sides of the located species. Note that there is no noncovalent interaction isosurface for one halogen atom in S2F1–DP–X₂ (X = F, Cl, Br), which refers to lack of such interaction in this case (and existing covalent bond). However, there are much more troughs in the negative area of electron density for the S2F1–DP–X₂ (X = F, Cl, Br) complexes in Fig. 6, which could be assigned to more attractive interactions in these models.

Conclusion

This theoretical study focused on the noncovalent interactions of SFs with halogens and halides (with the exception of iodine) for possible sensor and adsorption applications. The initial SFs were chosen from optimized C₆₀-fullerene and by the substitution of sulfur atoms and locating halogens and halide on logical distances above the surface of the fullerene. Energy calculations showed that all the interactions are thermodynamically favorable, more significantly in the gaseous phase than the solvent phase. Among the solvents, in the cyclohexane phase, the interaction energy showed that it was the most favorable solvent and benzene and chloroform are placed 2^nd and 3^rd, respectively. Moreover, the obtained data showed that the chalcogen bond is the most important noncovalent interaction in these systems (with the highest interaction energies), the halogen bond stands 2^nd in importance and the π-halogen bond does not exist in most of these complexes. For the S2F1–DP–X₂ (X = F, Cl, Br) complexes, the solvent phase appears to be more desirable media for interaction compared to the gaseous phase, which refers to the fact that the mentioned complexes are more polar than their related initial structures. Donor–acceptor interaction energies, obtained from NBO calculations, demonstrated the direct effect of S-doping on acceptor species as long as σ* C–S is anticipated as the strongest acceptor. However, more powerful transitions from the lone-pair of the carbon next to the doped sulfur to the σ* of halogens were observed. NBO atomic charge revealed that the doped sulfur atom has a large positive charge, especially when it interacts with halides, and partial positive charge for neighborhood carbon atoms in the doped site of fullerene. The investigation of DOS plots and population data shows a decrease in E_g upon S-doping of fullerene, chemical hardness tended to lower values, global softness was increased and the electrophilicity index decreased. The smallest energy gap was observed for the S2F1 fullerene. QTAIM information indicated that these noncovalent interactions would exist along the absorbent and adsorbate for complexes with electrostatic type interactions. Interestingly, the S2F1–DP–X₂ (X = F, Cl, Br) complexes have some other type of interactions due to their ρ values and negative Laplacian, which likely shows the partial functionalization of the S2F1 fullerene because of its particular doping and special locating of halogens above it. It has been proven that in these cases, the [S2F1–X]⁺X₃⁻ stable structure could be produced from S2F1 and two molecular halogens via an exothermic process. NCI index calculations confirmed the noncovalent interactions via depicted isosurfaces, and some troughs appeared in the sign λ₂ × ρ versus RDG plots, which are proof for the occurrence of electrostatic noncovalent interactions. Moreover, it is noteworthy that there is no noncovalent isosurface for one halogen atom that binds with the carbon atom next to the doped sulfur atom in the S2F2–X2 complexes.

Acknowledgements

We are grateful to High Performance Computing Research Center (HPCRC) at Amirkabir University of Technology, Tehran, Iran (National Amirkabir supercomputer; http://hpcrc.aut.ac.ir) and National High-Performance Computing Center (NHPCC) at Isfahan University of Technology (http://nhpcc.iut.ac.ir, Rakhsh supercomputer) for providing computational facilities for performing this work. This work has been supported by research affair of Isfahan University of technology.

Notes and references

1. H. W. Kroto, J. R. Heath, S. C. O'Brien, R. F. Curl and R. E. Smalley, Nature, 1985, 318, 162–163.
2. R. E. Smalley, Acc. Chem. Res., 1992, 25, 98–105.
3. M. S. Dresselhaus, G. Dresselhaus and P. C. Eklund, Science of Fullerenes and Carbon Nanotubes, Academic Press Ltd., 1996.
4. A. Mittelbach, W. Hönle, H. G. von Schnering, J. Carlsen, R. Janiak and H. Quast, Angew. Chem., Int. Ed. Engl., 1992, 31, 1640–1642.
5. R. C. Haddon, A. F. Hebard, M. J. Rosseinsky, D. W. Murphy, S. J. Duclos, K. B. Lyons, B. Miller, J. M. Rosamilia, R. M. Fleming, A. R. Kortan, S. H. Glarum, A. V Makhija, A. J. Muller, R. H. Eick, S. M. Zahurak, R. Tycko, G. Dabbagh and F. A. Thiel, Nature, 1991, 350, 320–322.
6. M. D. Diener and J. M. Alford, Nature, 1998, 393, 668–671.
7. G. P. Tegos, T. N. Demidova, D. Arcila-Lopez, H. Lee, T. Wharton, H. Gali and M. R. Hamblin, Chem. Biol., 2005, 12, 1127–1135.
8. A. Jorio, G. Dresselhaus and M. S. Dresselhaus, Carbon Nanotubes: Advanced Topics in the Synthesis, Structure, Properties and Applications, Springer, 2007.
9. M. S. Dresselhaus and G. Dresselhaus, Adv. Phys., 1981, 30, 139–326.
10. Z. Wang, M. Yao, S. Pan, M. Jin, B. Liu and H. Zhang, J. Phys. Chem. C, 2007, 111, 4473–4476.
11. S. Glenis, S. Cooke, X. Chen and M. M. Labes, Chem. Mater., 1996, 8, 123–127.
12. H. Jiao, Z. Chen, A. Hirsch and W. Thiel, Phys. Chem. Chem. Phys., 2002, 4, 4916–4920.
13. H. Tavakol and F. Hassani, Struct. Chem., 2015, 26, 151–158.
14. H. Tavakol and A. Mollaei-Renani, Struct. Chem., 2014, 25, 1659–1667.
15. P. A. Denis, R. Faccio and A. W. Mombru, ChemPhysChem, 2009, 10, 715–722.
16. F. Hassani and H. Tavakol, Sensors Actuators B Chem., 2014, 196, 624–630.
17. E. G. Lewars, Computational Chemistry: Introduction to the Theory and Applications of Molecular and Quantum Mechanics, Springer, 2nd edn, 2011.
18. L. M. Azofra S. Scheiner, J. Chem. Phys., 2015, 142, 034307.
19. X. Guo, X. An and Q. Li, J. Phys. Chem. A, 2015, 119, 3518–3527.
20. M. D. Esrafili and F. Mohammadian-Sabet, Chem. Phys. Lett., 2015, 628, 71–75.
21. H. Tavakol, J. Phys. Chem. A, 2013, 117, 6809–6816.
22. H. Tavakol, Int. J. Quantum Chem., 2011, 111, 3717–3724.
23. M. D. Ganji, N. Ahmadian, M. Goodarzi and H. A. Khorrami, J. Comput. Theor. Nanosci., 2011, 8, 1392–1399.
24. Y. L. Wang, K. H. Su and J. P. Zhang, in Advanced Materials Research, Trans Tech Publ, 2012, vol. 463, pp. 1488–1492.
25. Y. F. Zhang, Y. F. Wang, N. Chen, Y. Y. Wang, Y. Z. Zhang, Z. H. Zhou and L. M. Wei, Nano-Micro Lett., 2010, 2, 22–25.
26. H. Zeng, Z. Sun, Y. Segawa, F. Lin, S. Mao and Z. Xu, Chem. Phys. Lett., 2001, 338, 241–246.
27. H. Tavakol, J. Mol. Struct.: THEOCHEM, 2010, 954, 16–21.
28. H. Tavakol, M. Esfandyari, S. Taheri and A. Heydari, Spectrochim. Acta, Part A, 2011, 79, 574–582.
29. C. W. Bauschlicher Jr. and H. Partridge, Chem. Phys. Lett., 1995, 240, 533–540.
30. H. Tavakol, Struct. Chem., 2011, 22, 1165–1177.
31. H. Tavakol, J. Mol. Struct.: THEOCHEM, 2010, 956, 97–102.
32. H. Tavakol, Mol. Simul., 2010, 36, 391–402.
33. J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615–6620.
34. J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2008, 128, 084106.
35. J. Antony and S. Grimme, Phys. Chem. Chem. Phys., 2006, 8, 5287–5293.
36. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
37. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3094.
38. G. E. S. M. J. Frisch, G. W. Trucks and H. B. Schlegel, et al., GAUSSIAN 09 (Revision A.02), Gaussian Inc., Wallingford, CT, 2009.
39. H. Tavakol and D. Shahabi, J. Phys. Chem. C, 2015, 119, 6502–6510.
40. F. Weinhold and C. R. Landis, Chem. Educ. Res. Pract., 2001, 2, 91–104.
41. N. M. O'boyle, A. L. Tenderholt and K. M. Langner, J. Comput. Chem., 2008, 29, 839–845.
42. R. G. Parr, R. A. Donnelly, M. Levy and W. E. Palke, J. Chem. Phys., 1978, 68, 3801.
43. R. F. W. Bader, Atoms In Molecules: A Quantum Theory, Oxford University Press, 1990.
44. T. A. Keith, AIMAll, Version 10.05.04, TK Gristmill Software, Overland Park KS, USA, 2010.
45. E. R. Johnson, S. Keinan, P. Mori-Sanchez, J. Contreras-Garcia, A. J. Cohen and W. Yang, J. Am. Chem. Soc., 2010, 132, 6498–6506.
46. J. Contreras-García, E. R. Johnson, S. Keinan, R. Chaudret, J.-P. Piquemal, D. N. Beratan and W. Yang, J. Chem. Theory Comput., 2011, 7, 625–632.
47. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
48. T. Williams, K. Colin, et al., gnuplot, version 4.6.5, http://www.gnuplot.info, 2014.
49. J. Cerny and P. Hobza, Phys. Chem. Chem. Phys., 2007, 9, 5291–5303.
50. J. R. Lane, J. Contreras-García, J.-P. Piquemal, B. J. Miller and H. G. Kjaergaard, J. Chem. Theory Comput., 2013, 9, 3263–3266.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5ra08141f

---