Theoretical study of the interaction of fullerenes with the emerging contaminant carbamazepine for detection in aqueous environments

The global increase in drug consumption exposes the growing need to develop new systems for the detection, capture, and treatment of bioactive molecules. Carbamazepine is one instance of such contaminants at the top of the ranking commonly found in sewage treatment systems. This work, therefore, presents a theoretical study of fullerene C60 and its derivatives with substitutional doping with B, Al, Ga, Si, Ge, N and P, for the detection and capture of carbamazepine is aqueous medium. Solvation effects were included by means of the Polarizable Continuum Solvent method. The results indicate that doped fullerenes are sensitive for the detection of carbamazepine both in gaseous and aquatic environments. Investigation on the intermolecular interactions between the drug and the fullerene molecule were carried out, allowing the characterization of the interactions responsible for stabilizing the adsorption of carbamazepine to the fullerenes. The theoretical survey revealed that fullerenes doped with Al, Ga, Si and Ge chemically adsorb carbamazepine whereas for the case of fullerenes doped with other heteroatoms physisorption is responsible for the molecular recognition. Relying on DFT calculations, the fullerene derivatives C59Al, C59Si and C59Ga are the most suitable to act both as a sensor and to uptake carbamazepine in aquatic environments.


Methodology and computational details
Geometric and electronic properties of the doped-fullerenes, CBZ and the respective dimers were determined under the framework of Density Functional Theory (DFT) employing the ωB97XD exchange-correlation functional coupled with the 6-31G(d) basis set. Pople-type of double-ζ quality 6-31G(d) basis set was employed for it being a suitable choice to deal with relatively large systems with a good compromise between accuracy and computational cost [35][36][37][38] and it provides fairly good geometries and energetics when compared with experimental results, especially for large systems 24,35,36,[39][40][41] . Vibrational frequencies were calculated at the same level of theory to confirm that all structures correspond to genuine minima. These calculations were performed with the quantum chemistry Gaussian 09 42 suite of programs. Solvent effects were included using the Polarizable Continuum Model (PCM) 34 . Cartesian coordinates of the investigated species can be found in the Supplementary Information file, Tables S1, S2 and S3. The functional ωB97XD 43 is a long-range-corrected hybrid functional that includes Grimme´s empirical dispersion model D2 44 , which adds an extra term to the energy obtained by DFT calculations. ωB97XD is a good choice to include London dispersion interactions which are crucial for the dimers studied here and shows good performance to treat non-covalent interactions 38 . Besides that, our previously study showed that ωB97XD brings a good cost benefit for calculations involving molecular adsorption by fullerenes 37 . Concerning HOMO-LUMO gap energies (E g , HOMO-Highest Occupied Molecular Orbital and LUMO-Lowest Unoccupied Molecular Orbital), we selected M06L functional 45 , relying on previous theoretical investigations revealing that this functional is suitable to predict the experimental value of E g for the C 60 fullerene 18 .
Adsorption energies ( E ads ) were computed within the supramolecular approach, in which the energy of the complex ( E complex ) is subtracted from the energies of the monomers ( E CBZ + E C 59 X ), All energetics were corrected with the Basis Set Superposition Error (BSSE) using the counterpoise method 46 . Partial atomic charges were evaluated using ESP 47 charge analysis and from the molecular electrostatic potential (MEP) map of the systems were plotted.
To investigate the type of intermolecular interaction, Quantum Theory of Atom and Molecules (QTAIM) [48][49][50][51][52] and Reduced Density Gradient (RDG) index 53,54 were employed. Both approaches are based in the analysis of the electron density (ρ ) and the second eigenvalue ( 2 ) of the Hessian matrix (∇ 2 ρ = 1 + 2 + 3 ) , at the called Bond Critical Points (BCPs) in the region between interacting molecules. QTAIM and RDG analysis has been applied to help categorizing non-covalent interactions in different systems 33,[55][56][57] . QTAIM and RDG properties were computed using the wavefunction analysis free program Multiwfn 58 . The plot of the isosurfaces and molecules for both QTAIM and RDG analysis were made with VMD software version 1.9.3 59 .

Results and discussion
Doping effects on the frontier orbitals. To investigate the electronic properties of the C 59X fullerenes, we employed the M06L/6-31G(d)//wB97XD/6-31G(d) model. The M06L exchange-correlation functional is well-suited to predict HOMO-LUMO gap energies 18 in close agreement with the experimental E g values ranging from 1.80 60 up to 1.90 eV 61 in the case of pristine fullerene (see Table 1).
The C 60 Highest Occupied Molecular Orbital (HOMO) showed quintuple degeneration and the Lowest Occupied Molecular Orbital (LUMO) showed a triple degeneration. As already observed in previous studies 62 , C 60 doping tends to break this degeneration. The reduction of the energy gap between the border orbitals is recurrently observed in doped carbon nanostructures 19,21,24,[63][64][65][66] . The results obtained for the C 59 X confirm this trend. Glancing at Table 1 and Fig. 1, it can be observed that, in general, the C 59 X exhibited a decrease in the energy www.nature.com/scientificreports/ gap, E g , relative to the C 60 fullerene in both vacuum and water. Figure 1 also portrays the marginal impact of the implicit solvation treatment on the energies of the frontier states and the energy gaps. The value of E g can be adopted as an index of chemical stability, so that lower values of E g is commonly associated with higher reactivity 67 . Thus, the energetic variation of frontier orbitals indicates an increase in chemical sensitivity of heterofullerenes 63 . The decreasing order of E g can be observed in Fig. 1 and this trend was ubiquitous for all doped fullerenes regardless of environment. Moreover, it can be observed that the reduction of E g in C 59 X, is more strongly influenced by the increase of orbital HOMO energy, i.e., E H , destabilizes whereas E L , was Table 1. Energies (in eV) of the frontier orbitals and HOMO-LUMO gap for the doped fullerenes and its percentual variation upon complexation of CBZ ( E g in% ) as determined at the M06L/6-31G(d)//wB97XD/6-31G(d) level of theory.

M06L/6− 31G(d)
Vacuum Water    68 . E H increases in the following order: C 60 , C 59 Ge, C 59 Si, C 59 B, C 59 Ga, C 59 Al, C 59 P and C 59 N, indicating that the electronic perturbation of the system, either by of electron rich or deficient-elements in the C 60 valence shell, produces destabilization of the HOMO orbital. For the elements of group 13, the insertion of the heteroatoms generates an electronic vacancy in the HOMO orbital. Hence, the formation of a Singly Occupied Molecular Orbital (SOMO) with higher energy than the C 60 HOMO orbital occurs 62,69,70 . The new SOMO orbital is positioned energetically above the HOMO orbital of neat fullerene. In the case of N-and P-doping, which has one extra electron in the valence shell, also leads to the formation of a new SOMO orbital [69][70][71] , which is also energetically positioned above of the HOMO orbital of pristine fullerene. For the C 59 Si and C 59 Ge, since they are isoelectronic in the valence-shell to the carbon atom, HOMO was less affected by the doping element (see lower panel of Fig. 1). For these elements the E g reduction is attributed by the more important decrease in E L . Such effect was already reported in the literature 18,63,68,72 . MEP and ESP partial atomic charges. Figure 2A brings the molecular electrostatic potential (MEP) maps for the doped fullerenes. Table S4 shows the values of ESP charges of the heteroatoms of each C 59 X structure.
For doped fullerenes, with B, Al, Ga and Si, an electrophilic behavior was observed in the region around the doping atom, as can be seen in Fig. 2A. For the fullerenes doped with Ge, N and P, one can see a nucleophilic behavior over the dopant atoms with charge accumulation. It is also possible to note a positive charge region in the carbons near the Ge, N and P atoms. Through the MEPs for the doped fullerenes, it is observed that the heteroatom generates a significant non-local electronic perturbation over the fullerene cage. It is also noted that, around the dopant, there is an electrophilic region with can represent an important site for interaction with CBZ. Still looking at Fig. 2A, it is noted that the drug has a reactive region in the amide group and in the below the aromatic rings opposite to the amide moiety displaying nucleophilic behavior. In view of the predominantly electrophilic characteristics of doped fullerenes, it is assumed, therefore, that the CBZ molecule will have a greater tendency to interact with C 59 X along the carbonyl fragment.
Inspired by the MEPS displayed in Fig. 2A, we performed the optimization of the complexes CBZ + C 59 X considering two configurations of the CBZ molecule as displayed in Fig. 2B. The configuration (I) was considered, www.nature.com/scientificreports/ observing the largest superficial area of the aromatic fragments to establish interactions between the fullerenes and the drug. Thus, it is expected both configurations will display non-negligible interactions and hence, configurations (I) and (II) were conducted in this work.

Frontier orbitals energy variations upon complexation of CBZ.
The sensing capability of doped fullerenes for CBZ adsorption was studied by considering the variation of E g after the adsorption of CBZ. Figure 3 depicts the effect of doping and interaction with CBZ of the E g and the percentage variation of the gap, E g , which measures the percentage variation of E g of the C 59 X after the CBZ adsorption on the fullerene cage (Table S5 present the numerical values of E H , E L , E g and �E g ) . The change of the E g is an important parameter, since the energy gap is directly correlated with electrical conductivity (σ), σ ∝ exp(− E g /2k B T) 73 (T is the temperature, and k B represents Boltzmann constant) so that if the E g of the doped-fullerene is strongly affected by the presence of CBZ an electronic signal enables identification of the drug. Figure 3 (top panel) portrays interesting features, such as, (i) similar to that observed in the case of isolated C 59 X, dopants decreases the E g gap (recall Fig. 1), (ii) PCM solvation treatment exhibit little impact on E g to order of ~ 0.1 eV, and (iii) the relative orientation of CBZ with respect to the doped site on fullerene, the relative configurations I and II, display significant differences only when the doping atom is B, Si and to a lesser extent Ge. To assess more profoundly the changes of E g when C 59 X interacts with CBZ, we can see from Fig. 3 (lowest panel), for configuration II in water, that the doping atoms B, Ga, Al, Si and Ge provokes large variations of the E g gap, that is, from an electronic standpoint, these heterofullerenes show more sensibility towards CBZ. Still focusing on Fig. 3 (lowest panel), we note that dopants N, P and Ge, the variation of the E g gap after complexation in vacuum, was positive, that is the interaction actually increased the gap afterwards complexation especially in configuration II (positive values of E g ). These effects might be attributed to electrostatic repulsions between the polar amide fragment of CBZ and the electron-rich dopants in C 59 X as discussed in the MEPs presented in Fig. 2.
Adsorption energies in gas and aqueous phase. Because ωB97XD is reported to deliver good interaction energies between system where non-covalent interactions are important 37,38,72,74,75 , we will henceforth discuss the interaction energies between CBZ and C 59 X relying on calculations performed at the ωB97XD/6-31G(d) level of theory. In Table 2 and Fig. S2, the results for the adsorption energy between the C 59 X molecules and the CBZ are presented for both vacuum and aqueous environments and were BSSE-corrected. The pristine C 60 showed the smallest adsorption energy in both configurations and environments. In general, the doping effect was to increase (in modulus) the E ads energy. Solvation effects had little impact on E ads , but for most complexes the interaction was enhanced in water environment. www.nature.com/scientificreports/ As can be observed in the optimized systems presented in Fig. S1 and in previous studies 21,24,62,63 , C 60 doping generates a structural deformation localized in the region of the heteroatom. Thus, it is plausible to assume that the atomic radius of the X atom influences the order of interaction between the C 59 X and the CBZ. An estimative for the atomic radius, of the heteroatoms, can be obtained through the medium variation of the bonds between the neighboring carbons to the dopant atom.
Due to the larger atomic radius, the atoms of Al, Ga, Ge, Si, P and B are projected outside of the C 59 X structure, which makes these heteroatoms more susceptible to nucleophilic/electrophilic attacks. This heteroatom projection can be observed by the variation of the bond length between the carbon atoms and the X atom. Table S4 presents the atomic radius 76 , R α , of the heteroatoms in the doped fullerene and the average bond length, R ave , between the heteroatoms and the surrounding carbon atoms. The R ave values obtained here agree with values obtained in previous studies 63,68 . Through the data presented in Table S4, it can be noted that the order of interaction, tends to follow the increasing order of R α and R ave . However, it is important to point out that, the electronic behavior of the C 59 X also may influence the interaction with CBZ.
Although the systems doped with P and Ge presented larger R α and R ave , the nucleophilic behavior of the C 59 Ge and C 59 P, at the doping site, may explain the higher value of E ads obtained for others C 59 X-CBZ in comparation with the values of E ads obtained for C 59 Ge-CBZ and C 59 P-CBZ. Besides that, it is also possible to note that, although Al and Ga present very similar R α and R ave values, Al atom has a higher ESP charge than Ga atom, as presented in Table S4. For this reason, the E ads is larger for C 59 Al than C 59 Ga. For both configurations, C 59 X interact with nucleophilic regions of CBZ. Therefore, it is expected that systems with more electrophilic characteristics, around the X atom, and with higher R α and R ave values, have a better interaction with CBZ. Thus, because of the more electrophilic behavior, presented in Fig. 1 and Table S4, and due the large value of R α and R ave , the C 59 Al system presents a higher interaction with CBZ in both, gas and aqueous environment.
However, an interesting exception to this tendency can be observed for the E ads obtained by C 59 Ga and C 59 Si, in conf. II and in both phases. Although Ga presented higher ESP charge, R α and R ave values, C 59 Si presents higher E ads than C 59 Ga. In order to understand this exception, the charge distribution, along the hexanes and pentane face surrounding the heteroatoms, was evaluated. It is well known that dominant partial charges, either positive or negative, accumulated on specific parts of interacting molecules have important influence for effective interaction 74,77,78 . Bearing it in mind, the sum of the ESP charges in the two hexanes faces and in the pentane face around the doping atoms was taken into account (See Fig. 2C). These results are presented in Table S4. By means of these data, it was observed that, the C 59 Si has the highest positive value of the total charge in the pentane face. For this reason, CBZ reoriented to interact specially with this face. This pentane-CBZ interaction is responsible for increasing of E ads . This interaction, together with a possible electrostatic dipolar bond, discussed in more details in "Non-covalent interaction" section, would be enough to stabilize C 59 Si when compared to C 59 Ga.
To clarify the nature of these interactions in the C 59 X-CBZ systems, a detailed study of the intermolecular interactions, in both configurations and phases, was carried out.

Non-covalent interactions.
To investigate the nature of the interaction between doped fullerenes and the CBZ drug, the Reduced Density Gradient (RDG), together with Quantum Theory of Atoms in Molecules www.nature.com/scientificreports/ (QTAIM) were employed. In the RDG analysis, the combination of electron density ρ and the signal of 2 provides a tool to distinguish between repulsive or attractive interactions. For this task, the function sign ( Tables S8 and S9, are represented by orange points and the bond paths by the yellow lines. By means of the data presented in Tables S6-S9, one can see that, according to the QTAIM scheme, in configuration (I) almost all systems presented van der Waals interaction. The exceptions were observed for Al and Ga which presented interactions with character resembling polar covalent bonds, in both gas and aqueous phase.
For C 59 Al-CBZ (II), in both phases, the values observed in Tables S6 and S7 indicate that an interaction with ionic character was observed between Al and O from CBZ. For C 59 B-CBZ (II), C 59 Ga-CBZ (II), C 59 Si-CBZ (II) and C 59 Ge-CBZ (II), in both phases, a partial covalent character interaction was observed between the heteroatom and the oxygen in CBZ. For the C 60 -CBZ (II) C 59 P-CBZ (II) and C 59 N-CBZ (II), for both phases, QTAIM parameters indicate that the intermolecular interactions for these complexes are exclusively of van der Waals type.
As mentioned before, the N and P dopants, in C 59 N and C 59 P, have a nucleophilic character which induces a reorientation of CBZ in front of the fullerene cage. Tables S8 and S9 present the optimized structure, in gas phase and aqueous environment respectively, where it can be seen this CBZ reorientation. For this reason, the CBZ in C 59 N-CBZ and C 59 P-CBZ, prefers interacting to the fullerene cage instead to the dopant atoms. The CBZ reorientation was observed in both phases and configurations. This may explain the exclusivity of van der Waals interaction in these complexes.
Additionally, in gas phase, for C 59 Si and C 59 Ge, one C · · · H electrostatic dipolar bonding was observed in each one of these systems. As discussed in the "Adsorption energies in gas and aqueous phase" section, this electrostatic dipolar C · · · H bonding (See values of BCP 1 for Conf. II at Tables S6 and S7 for C 59 Si-CBZ) together to the ESP accumulation charge, would help to explain the increase in the stabilization for C 59 Si-CBZ, comparing with the stabilization energy of C 59 Ga-CBZ (See Fig. 2 and Table 2). For C 59 Ge-CBZ, the value of E ads is affected by the nucleophilic character of Ge atom in the fullerene cage (See Table S4). For this reason, the C · · · H electrostatic dipolar bonding does not produce a notable increase in the stabilization of C 59 Ge-CBZ system. Now focusing on RDG analyses, the isosurfaces were associated with the scattering graphs, shown in Fig. 4 and Tables S8 and S9. In Fig. 4 the RDG isosurfaces and the scatter graph for C 59 Al-CBZ is presented. Greencolored isosurfaces indicate non-covalent interactions while red-ones is due to steric effects. The blue-colored www.nature.com/scientificreports/ regions represent strong attractive effects. The RDG isosurfaces and scatter graph for all C 59 X-CBZ complexes are detailed in Tables S8 and S9. For the C 59 B-CBZ and C 59 N-CBZ complexes, in both phases and configurations, the RDGs isosurfaces and scatter graphs are similar to the RDGs surfaces and scatter graphs observed for C 60 -CBZ. This may be due the heteroatom radius. B and N have R α close to carbon, the fullerene cage is less deformed and the contact area for vdW interactions is similar to the case of pristine fullerene. On another hand, it was observed that the scatter graph presented more points in the attractive region for C 59 Al-CBZ (I), C 59 Si-CBZ (I), C 59 Ga-CBZ (I) and C 59 Ge-CBZ (I), in both phases. Again, it is observed that for those larger heteroatoms NCI increases accordingly. These observations are in line with the interaction strength order of the C 59 X-CBZ complexes.
For these systems in conf. (II), QTAIM results indicated that the heteroatom-oxygen interaction may have a greater importance for E ads . However, as showed for both, QTAIM and RDG, NCIs still play a significant role in the stabilization of the systems, especially for C 59 Al-CBZ (II), which has an important vdW interaction in the green region of the RDG isosurfaces (see BCPs 2 in Tables S6, S7, S8, S9 and Fig. 4), and for C 59 Si-CBZ and C 59 Ge-CBZ (II), which have an indication of an electrostatic interaction (see BCPs 1 in Tables S6 and S8).
By means of non-covalent interaction analyses, it was observed that, for all complexes in configuration (I), and for C 60 -CBZ (II), C 59 P-CBZ (II) and C 59 N-CBZ (II), the drug is physisorbed to the fullerene cage. Nevertheless, for C 59 Al-CBZ (II), C 59 B-CBZ (II), C 59 Ga-CBZ (II), C 59 Ge-CBZ (II) and C 59 Si-CBZ (II), a chemical interaction occurs between CBZ and these doped fullerenes. We should note, however, that C 59 Ga-CBZ (II) exhibited another local minimum, and this configuration is presented in Fig. 5.
Aiming to find the true minimal structure for C 59 Ga-CBZ (II) a relaxed scan was made having as starting point the structure presented at the left of Fig. 5. For this relaxed scan, we rotated the CBZ in front of the C 59 Ga, making a variation of 10° for each position. All 36 calculations converged to the structure presented at the right in Fig. 5 which indicated that this is the true minimal configuration for C 59 Ga-CBZ (II).
One can see that, comparing this local minimum region to the global minimum region, the presence of an extra vdW interaction helped in the increase of E ads in 0.016 (eV). This small increase in the energy confirms that the interaction between the Ga atom to the O atom from CBZ may represent the major influence for the E ads in the C 59 Ga-CBZ (II). However, besides these new stable configurations, it is important to highlight that, due the small energy difference in the E ads , these possible new configurations, for C 59 X-CBZ (II), may not have influence on the CBZ adsorption capacity presented by the C 59 X molecules. The results indicate that, due their E ads values and electrophilic character, C 59 Al, C 59 Si and C 59 Ga are the doped fullerenes more indicated for acting as both, sensor, and capture system of carbamazepine. On the left, a local minimal structure presents the lack of a vdW interaction. On the right, after a new scan calculation, the missed vdW interaction appears. This new vdW interactions helps in a slight increase of the E ads . www.nature.com/scientificreports/

Conclusions
The main purpose of this study is to investigate the adsorption of carbamazepine in pristine and doped fullerenes to disclose its capability to sense this drug. Two relative orientations of CBZ with respect to the doping site were considered. Frontier orbital analysis showed that the HOMO-LUMO energy gap decreases whether employing heteroatoms from group 13, 14 or 15. Solvation effects modelled by means of the PCM approach had little impact on electronic, geometric, and energetic properties. Upon complexation, it was observed a gap variation in almost all C 59 X-CBZ systems. In water, almost all complexes showed a gap reduction after the CBZ interaction.
The results indicated that this gap variation would produce a signal due to changes in electrical conductivity of doped fullerene when CBZ adsorbs. The adsorption energies were enhanced for every heteroatom, but the most stable complex was obtained with Al-doping. Group 15 elements do not stabilize the complexes significantly and therefore, are not suitable to detect organic structures with polar functional groups such as the amide fragment of CBZ. Intermolecular interaction analyses were conducted within the framework of the QTAIM descriptors. It was observed that a chemical interaction occurs for C 59 Al-CBZ (II), C 59 B-CBZ (II), C 59 Ga-CBZ (II), C 59 Ge-CBZ (II) and C 59 Si-CBZ (II). Almost all these complexes presented a partial covalent character interaction between the doping atom and the oxygen in CBZ molecule. The exception was observed for C 59 Al-CBZ (II), which presented an ionic character interaction between the Al atom, from C 59 Al, and the O atom, from CBZ. For all complexes in configuration (I), and for C60-CBZ (II), C59P-CBZ (II) and C59N-CBZ (II), the drug has a physic adsorption. The cage structural and electronic deformation was evaluated to understand the interaction order. It was observed that, for those doped fullerenes with higher local structural deformation and higher electrophilic character, the drug adsorption was favored. In another word, the interaction order follows the heteroatom size and the electrophilic behavior around the doping region. This was corroborated by the analysis of NCI. Our calculations showed that C 59 Al-CBZ (II) presented the highest value for E ads with, − 2.411 eV in gas phase and − 2.697 eV in aqueous media. Finally, our results indicate that, due their greatest adsorption energy, in both phases and configurations, C 59 Al, C 59 Ga and C 59 Si are potential candidates aiming at drug sensor applications.

Data availability
The data presented in this study are available in the article and in the Supporting Information file. www.nature.com/scientificreports/