Non-covalent interactions of cysteine onto C60, C59Si, and C59Ge: a DFT study

The study of intermolecular interactions is of great importance. This study attempted to quantitatively examine the interactions between cysteine (C3H7NO2S) and fullerene nanocages, C60, in vacuum. As the frequent introduction of elements as impurities into the structure of nanomaterials can increase the intensity of intermolecular interactions, nanocages doped with silicon and germanium have also been studied as adsorbents, C59Si and C59Ge. Quantum mechanical studies of such systems are possible in the density functional theory (DFT) framework. For this purpose, various functionals, such as B3LYP-D3, ωB97XD, and M062X, have been used. One of the most suitable basis functionals for the systems studied in this research is 6-311G (d), which has been used in both optimization calculations and calculations related to wave function analyses. The main part of this work is the study of various analyses that reveal the nature of the intermolecular interactions between the two components introduced above. The results of conceptual DFT, natural bond orbital, non-covalent interactions, and quantum theory of atoms in molecules were consistent and in favor of physical adsorption in all systems. Germanium had more adsorption energy than other dopants. The HOMO–LUMO energy gaps were as follows: C60: 5.996, C59Si: 5.309, and C59Ge: 5.188 eV at B3LYP-D3/6–311 G (d) model chemistry. The sensitivity of the adsorption increased when an amino acid molecule interacted with doped C60, and this capability could be used to design nanocarrier to carry cysteine amino acid.


Introduction
The first attempts to come up with the concept of intermolecular interactions were made by Clausius [1], van der Waals [2], and London [3,4]. Intermolecular interactions can be divided into long-range and short-range classes. The longrange class includes electrostatic, induction, and dispersion forces, and it varies with the inverse powers of the distance rn , which is the reciprocal of the intermolecular distance. Conversely, the short-range class includes exchange and repulsion forces that decrease exponentially with distance, as in e -αr . Both repulsive and attractive electrostatic forces arise from classical Coulombic interactions between the charge distributions of two molecules. These forces are also pairwise additive and anisotropic. The interacting molecules cause instantaneous fluctuations in the electron distribution, and such a disruption creates dispersion forces that are pairwise additive and always attractive. Induction forces, which are non-additive and attractive, are created by the distortion of the distribution of molecular charges resulting from the influence of the electric fields of other molecules. Exchange and repulsion are both non-additive and of opposite signs [5].
Ab initio calculations are a reliable method for obtaining polarizabilities (i.e., a set of constants to show the charge redistribution when a molecule is exposed to an electric field) and electric multipole moments used to explain the long-range forces [6]. In the case of short-range forces, the theory is rather more complicated [7] due to the overlap of electron densities between molecules. The previous longrange theory should be modified in terms of short-range penetration [8][9][10], charge transfer [11], and damping effects [12][13][14]. Generally, full quantum methods can only be well implemented in small systems, but there are technical problems with large molecules [15]. The development of wave function analysis methods has helped predict the nature of intermolecular interactions, even in the case of large systems. Although the ab initio methods are few, neat, and classifiable, there are countless wave function analysis methods that have not yet been systematically categorized in the scientific literature.
This study examined the adsorption of cysteine amino acid molecule onto the surface of pristine-, silicon-, and germanium-doped fullerene nanocages, C 60 , C 59 Si, and C 59 Ge, respectively. All structures, such as isolated species and complex systems, were optimized using various functionals, including B3LYP-D3, ωB97XD, M062X, and 6-311G (d) basis set. The most stable structure in terms of electronic energy was chosen for further wave function analyses using the B3LYP-D3/6-311G(d) model chemistry. We referred to the conceptual density functional theory (DFT), natural bond orbital (NBO), non-covalent interactions (NCI), and quantum theory of atoms in molecules (QTAIM) from among the most important and reliable wave function analyses ever developed. In this study, they were used to understand the nature of intermolecular interactions. All studies were performed in vacuum; therefore, the sensitivity and reactivity of the adsorbent were appreciated. The main ideas of each analysis were explained in the shortest possible way, as these materials are scattered in the available resources, and it is a little difficult for beginners to have them together.

Computational details
The Kohn-Sham DFT framework was implemented in this study. Each structure, including amino acidand nanocages, was geometrically optimized in a vacuum using the so-called Berny geometry optimization algorithm developed by H. B. Schlegel in 1982 [47]. The PBE0 functional [48,49], which is the hybrid-exchange correlation form of the Perdew-Burke-Ernzerhof generalized gradient approximation (GGA) [50], meta-hybrid GGA functional M06-2X developed by the Truhlar group [51,52], the Head-Gordon group functional ωB97XD [53], which includes dispersion and long-range corrections, and the Becke 3-parameter Lee-Yang-Parr functional with Grimme dispersion correction B3LYP-D3 [54][55][56] were employed. The MC-311 or 6-311G (split-valence triple-zeta) basis set in the Gaussian package [57] has been developed by many individuals [58][59][60][61][62][63][64][65][66]. In the current work, the 6-311(d) basis set was used (i.e., d-type Cartesian-Gaussian polarization functions). Split-valence basis sets allow orbitals to alter the size but not the shape. When polarized functions are added to the basis sets, this limitation is removed by including orbitals with an angular momentum greater than what is needed for the ground state to depict each atom. This basis set is large enough and well able to simulate molecular orbitals. Data from benchmark studies also confirm this fact [67][68][69][70][71]. GaussView 6.0.16 [72] and ChemCraft [73] packages were used for building molecules. Linux-based Gaussian 16 software Rev. C.01 [57] was used for self-consistent field (SCF) calculations (as link 502). The default convergence criteria remained intact during the calculations. The SCF convergence was considered by comparing the maximum force and maximum displacement with the threshold values (0.00045 Hartree/Bohr and 0.0018 Bohr). Moreover, no symmetry limitations were imposed on the optimization process. Wave function stability, frequency checks, and zero-point energy corrections (ZPECs) were considered to ensure the accuracy of the calculations. Stability calculations guarantee that this optimized electronic wave function may be minimal in the wave function space rather than a saddle point and that it is completely different from finding minima or saddle points on a nuclear potential energy surface. As the Gaussian software includes the NBO version 3.1 software [74][75][76] (as link 607), it is used to perform population analysis studies. The Multiwfn [77] package developed by Tian Lu is fed by formatted Gaussian checkpoint files for various wave function analyses. In this work, Multiwfn was used for NBO, NCI, and QTAIM studies. O'Boyle et al. developed the cclib and GaussSum [78] package to obtain the DOS diagrams.
To calculate the adsorption energy (E ads ) of the two molecules (nanocage and cysteine), the following underlying relation is applied: where E cage/amino is the energy of the cluster, and E amino and E cage are the energies of the isolated amino acidand nanocage, respectively. The negative values of E ads (i.e., exothermic adsorptions) show that the formed amino acid/ nanocage cluster is stable. The electron density of each nucleus can be determined using a function centered on another nucleus. Therefore, in all structures, the quality of (1) the basis set is not the same. This means that the basis set of one molecule can be effective in compensating for the violation of the basis set of another molecule. This effect is called the basis set superposition error (BSSE). There are two well-known methods for the BSSE correction, (1) chemical Hamiltonian approach (CHA) [79] and Boys and Bernardi's counterpoise (CP) correction procedure [15,80,81] which CP method is employed in this work: The ZPECs are calculated using Eq. 3:

Geometric surveys
In the first step, molecular geometries were constructed using GussView software. The method of making nanocages is explained fully in the Help section of the software. Creating a cysteine amino acid molecule does not take much effort, but for more certainty, geometry related to amino acidmolecules was prepared using the PubChem online database [82]. At the beginning of the research, we chose models that were faster and less expensive. All structures, including isolated molecules and clusters, were initially examined using PBE0/6-311G(d) model. As our ultimate goal was to examine intermolecular forces, we used models that calculated intermolecular interactions more accurately. The optimization process was repeated with the same basis set and functionals M06-2X, B3LYP-D3, and ωB97XD. The chemistry of the different models varied within the trade-offs made between computational cost and accuracy.
Once geometric optimization was completed, the presence or absence of imaginary frequencies (i.e., frequencies that have negative values do not represent the minimum and actually represent a transition state) was examined. Thus, in each optimized structure, frequency calculations were performed. In the output file of the Gaussian software, the keyword Nimag was assigned for this purpose. A value of zero (Nimag = 0) means that there are no negative frequencies according to the number of negative eigenvalues of the Hessian matrix. Frequency calculations were performed to determine the ZPEC, which must be added to the total energy.
To study the interactions between nanocages and an amino acid molecule, a C 60 cage with 1.45-Å bond length was selected. The optimized C60 structure was used to prepare a doped nanocage. For this purpose, doped nanocages were obtained by substituting silicon and germanium instead of carbon atom and they were re-optimized in the same optimization process.
Each nanocage is illustrated in Fig. 1, which shows the position of the doped elements compared with the pure nanocage. The injection of the doped elements slightly changed the length of the bonds, indicating that the electronic structure was altered. These changes in the electronic structure led to differences in properties. Doping is a critical and successful strategy for detecting the properties of a nanomaterial as an adsorbent [83]. The dopants alter the sensing properties by changing the HOMO-LUMO gap and morphology, creating more centers for amino acid interaction on the adsorbent surface.
As shown in Fig. 2, there are five sites on the nanocage that have the potential to absorb the amino acid molecule. The top of the carbon atom is represented by T 1 . The position of T 2 is related to the placement of the amino acid molecule on top of the bond between two hexagonal rings. The T 3 position is the placement of the amino Fig. 1 The values of bond length for a C 60 , b C 59 Si, and c C 59 Ge. The optimization process has been done using the B3LYP-D3/6-311G (d) level of theory acid molecule on top of the bond between hexagonal and pentagonal rings. The T 4 position is located above the hexagonal ring. And T 5 is on top of the pentagonal ring. Structurally, the cysteine molecule has four heads, each of which can be placed on any of the five adsorption sites in C 60 (Fig. 3).
The starting point of the study is placing the amino acid molecule on the absorbent surface in different positions with different angles and distances and with different orientations and optimizing the whole cluster system. Finding the global minimum is the most challenging step of the study. Although we may never get it right away in the computing process, various strategies can be used to find many local minima and choose the most stable mode based on their quantitative value. If we have the full potential energy surface (PES) of two fragments, we can easily report the global minimum. However, having the PES requires thousands of calculations, which are practically impossible. The approach followed in this research is based on placing the amino acid molecule on the adsorbent surface at different distances and orientations. For this purpose, a low-cost method, such as PM6, was used. Dozens of initial orientations were selected to achieve the most stable structure in terms of energy content. The PEB0/6-311G(d) model chemistry was used to optimize the obtained stable structure and repeat the process using the abovementioned functionals.
The relaxed cluster structures obtained from B3LYP-D3/6-311G (d) are illustrated in Fig. 4. When the value of the energy of adsorption is below the range of chemical interest (i.e., from the third decimal place onwards following the decimal point), the results are identical [84]. Table 1 lists the values obtained from the calculations of the four different methods for optimizing complex structures. The quantitative amounts of adsorption energies obtained from the ωB97XD and B3LYP-D3 models were similar. Moreover, PBE0 data showed how different the data would be if the dispersion effect was not considered. Among the four different models, the values obtained from the B3LYP-D3/6-311G (d) model showed the highest stability, based on the absolute value of the highest absorption energy. As indicated in Table 1, among the doped elements, germanium produced the highest absorption energy (− 1.121 eV).
To understand the nature of intermolecular forces, whether they are strong or weak or are the result of van der Waals or electrostatic effects, the wave function is analyzed in the following sections. According to Table 1 and based on our previous experience [40,43], the B3LYP-D3 functional is a better match to the systems studied in this paper. Thus, it was used for the wave function analysis calculations.
Parr et al. [88] proved that μ is a quantitative measure of the tendency of an electron to pull out a system. The direction is from high to low μ values. The electronic chemical potential (μ) and electronegativity (χ) relation are expressed as Eq. 4 Pearson [89] showed hardness (η) as a measure of resistance of a system to a change in the electronic cloud as follows: Electrophilicity (ω) [90] is not the coefficient of the above-mentioned expansion. At the equilibrium point (μ = 0), ω is the stabilization energy gained by a system (i.e., ω is the capacity to accept an arbitrary number of electrons).
As shown in the above equations, the HOMO-LUMO energy gap (HLG) [91] is related to ionization potential (IP) and electron affinity (EA). Theoretically, based on the Koopmans' [92] and Janak's [93] approximations, the ionization potential is equal to the negative value of HOMO, HOMO = −IP , and the electron affinity is equal to the negative value of LUMO, LUMO = −EA . It should be noted that HOMO and LUMO can only be obtained from HF or DFT (i.e., single determinant methods) and the LUMO has no contribution to the total energy of the system. In addition, calculating the LUMO values is basis set sensitive; however, methods have been developed that can eliminate the LUMO dependence on the basis set [94].
The values of these descriptors are presented in Table 2. The energy gap (E g ) of C60 was 5.996 eV using the B3LYP-D3/6-311G (d) model chemistry, and the adsorption of cysteine on it reduced the energy gap to 5.067 eV.
The silicon-and germanium-doped nanocages also reduced the HLG values. Molecules with large HLG had high values of η, whereas those with low HLG had low η. Table 2 shows that the LUMO values became more negative after the adsorption of cysteine onto the nanocages; therefore, HLG was reduced. Reactivity of species was obvious due to the values of the descriptors η and μ. Conductivity arose from a high level of reactivity, and the charge transfer between the cages and adsorbate could be implemented through an chemical device as a carrier. To obtain a visual image of how the energy levels, especially HOMO and LUMO, are positioned, the total density of state maps is used in Fig. 5.

NBO analysis
Canonical molecular orbitals (MO) lead to delocalized descriptions of electrons derived directly from SCF methods [95]. Delocalization means that the MO obtained from several electrons belongs to adjacent atoms. As this is physically difficult to imagine, to create physical intuition, we had to use a method that would give us an understanding. As the wave function is invariant under unitary transformations, localization is possible. The NBO [74] method, Fig. 3 The values of bond length for cysteine molecule. The optimization process has been done using the B3LYP-D3/6-311G (d) level of theory In this study, NBO methods were utilized to measure the bond order. The most popular methods used to find the bond order are the Mulliken bond order analysis [96], the Mayer bond order [97][98][99], and the Wiberg bond index (WBI) in Löwdin orthogonalized basis [100,101]. Compared with Mulliken and Mayer's bond orders, WBI has less basis set dependence, especially the basis set that includes diffuse functions, and it provides more accurate results. The WBI values are reported in Table 3. Accordingly, we can conclude that C 59 Si and C 59 Ge adsorbents are more active materials in this study in adsorbing cysteine compared with pristine C 60 . The WBI values show that the interaction of the amino acid molecule with C 60 can be classified as a weak interaction. Conversely, the interactions between silicon-and germanium-doped C 60 with the cysteine molecule are stronger than the van der Waals interactions. The results of the WBI are in good agreement with the adsorption energies reported in Table 1.

QTAIM analysis
The idea that a molecule is a collection of atoms linked by a network of bonds is not directly related to quantum mechanics and is beyond theoretical definition. Bader's topological quantum theory of atoms in molecules (QTAIM) analysis [102,103] was developed to answer the question, "What is an atom in a molecule and how does one predict its properties?" Based on the topology of the electron density, the structure of a molecule is revealed by the stationary points of the electron density and the gradient paths that originate and terminate at these points.
The amount of ρ(r) in BCP determines the bond strength and the bond order. In the case of the covalent bond, the Laplacian is less than zero (∇ 2 ρ(r) < 0), and the values of ρ(r) are large, indicating that the charges are concentrated between two nuclei. The positive values of the Laplacian (∇ 2 ρ(r) > 0) and the low values of ρ(r) indicate that the charge dissipates in the distance between the two nuclei and that the interactions can be classified as a closed-shell type, which is related to ionic bonds, hydrogen bonds, and van der Waals bonds [104]. In particular, in the case of hydrogen bonding, if the electron density is in the range of 0.035-0.002 and if the Laplacian of electron density is in the range of 0.0139-0.002, the bond can be a hydrogen type.
The values of the Lagrangian kinetic energy G(r), potential energy density V(r), and energy density H(r) = G(r) + V(r) can also be helpful in identifying the type of interactions. For a covalent bond, H(r) < 0; H(r)/ρ(r) > > 0; and G(r) /|V|(r) < 0.5. For non-covalent interactions, H(r) > 0; H(r)/ρ(r) > 0; and G(r) /|V|(r) > 1. The virial theorem [105] suggests a relationship between G(r), V(r), and ∇ 2 ρ(r), Another descriptor is the bond elliptical index (ɛ), which is a good criterion for detecting conjugation and The interaction energy (E ads ) for C 60 , C 59 Si, and C 59 Ge with cysteine molecule from different terminals. All values are in (eV)  hyperconjugation. When ɛ is large, we have an elliptical structure that indicates that the π character is large. If ε = 0, the bond is cylindrical (double or triple) and very stable. Here, ε can be used to refer to the stability of the interactions and is defined as follows [106]: Considering the values shown in Table 4, the Laplacian electron density ∇ 2 ρ(r) is positive for all interactions between the cysteine amino acid molecule and all nanocages. Therefore, the interaction between amino acidand nanocages is detected as non-covalent. Figure 6 shows the BCP position. The presence of critical points between the amino acid molecule and the nanocages emphasizes the strong interactions between the two components. Among the nanocages, C 59 Si and C 59 Ge are more capable of adsorbing the amino acid molecule, as the values of G(r)/|V|(r) for both of them are between 0.5 and 1. This means that the interaction tends to be strong in van der Waals. The other descriptors, H(r) and H(r)/ρ(r), agree with these results. The values obtained from ε show stable intermolecular interactions.

NCI analysis
The results obtained from the previous sections show that the interactions between the amino acid molecule and pure and doped nanocages are non-covalent. Therefore, to verify this, these interactions were examined in terms of NCI analysis to determine the accuracy of the results. NCI analysis uses two functions: (1) signλ 2 (r)ρ(r) (product of electron  Table 4 The AIM topological parameters, including electron density (ρ(r)), Laplacian of electron density (∇ 2 ρ(r)), the kinetic electron density G(r), potential electron density V(r), eigenvalues of Hes-sian matrix (λ), and bond ellipticity index (ε) at BCPs of the cysteine molecule and C 60 , C 59  density and the sign of the second Hessian eigenvector) and (2) reduced density gradient (RDG), which is a dimensionless form of a gradient of electron density [107].  Depending on the position of the second function in the diagram, three areas are created, indicating the type of interactions. In the signλ 2 (r)ρ(r) < 0 region, strong non-covalent interactions are found; in the signλ 2 (r)ρ(r) ≈ 0 region, relatively weak van der Waals interactions are defined; in the signλ 2 (r)ρ(r) > 0 region, repulsion forces are dominant [107,108]. Figure 7 compares the NCI plots for both isolated nanocages and amino acid/nanocage clusters. Repetitive results from previous analyses are also shown in the NCI analysis. That is, in the case of C 59 Si and C 59 Ge adsorbents, the adsorption intensity is stronger than in others. To reference this, sign λ 2 (r)ρ(r) ≈ 0 and RDG ≈ 0.5 should be considered.

Conclusion
Various wave function analyses were performed to study the intermolecular interactions between the cysteine amino acid and fullerene nanosorbents. The standard model chemistry in this study was B3LYP-D3/6-311G(d), which was implemented for the geometry optimization process and NBO calculations. The PBE0, ωB97XD, and M06-2X functionals and the 6-311G(d) basis set were also used for geometry optimization. Considering the dispersion in these functionals was a factor that completely changed the results of these calculations. In the optimization process, the spatial orientations of the two monomers relative to each other were the determining factors that led to the finding of the local minima. Spanning the entire potential energy surface for such interactions could not be attained. Therefore, it is important to use an appropriate algorithm to find the local minima. The comparison of the adsorption energy of different clusters showed that the adsorbents C 59 Si and C 59 Ge trapped the amino acid molecule more intensely. The results of QTAIM and NCI analysis identified the intermolecular interactions of the type of strong van der Waals interaction for these nanocages. As the amino acid and the mentioned adsorbents interacted well, these nanomaterials could be used in the design a nanocarrier for cysteine molecule.