Houdhaifa R Belhouchet

Tahar Abbaz

Amel Bendjedou

Abdelkrim Gouasmia

Didier Villemin

A computational study of the inclusion of β-cyclodextrin and nicotinic acid: DFT, DFT-D, NPA, NBO, QTAIM, and NCI-RDG studies

Keywords: Cyclodextrin, Vitamin B3, Density function theory, Inclusion complex

Forming complexes with cyclodextrins can protect nicotinic acid (vitamin B3) from premature metabolism and enhance the solubility and stability of this drug. In this work, the formation of the inclusion complex of the neutral form of nicotinic acid and β-cyclodextrin was achieved. The complex is modeled using PM3, PM6-D4H3, and PM7, by considering two orientations of the guest: A and B, one is from wide to narrow rim, and the second is from narrow to the wide rim, respectively. The global minima positions were re-optimized using three density function methods: MN-15, B3LYP, and PW6B95-D3 with polarized Pople basis set 6-31G(d) in gas and aqueous phase. Orientation A showed the minimum complexation energy where the carboxylic functional group of nicotinic acid is located on the primary hydroxyl rim of β-cyclodextrin and the pyridine ring is totally embedded in the cavity. To further our study on the nature of complexation and the interactions of this host-guest system, different calculations were done. The reactivity indices showed that orientation A is harder than B and more electrophilic; the charge transfer occurred from the host to the guest and was confirmed by the natural population analysis (NPA). The natural bond orbitals (NBO) reveal the delocalization of orbitals between the host and the guest, quantum theory of atoms in molecules (QTAIM) analysis, and non-covalent interaction (NCI) analysis based on a reduced density gradient (RDG) give a detailed description of the nature of interactions between the host and the guest such as the hydrogen bonding and van der Waals interaction, and confirmed the stability of the complex given by the orientation A.

Introduction

The host-guest complexes formed by the cyclodextrins have important uses in different fields: food industry, pharmacy, cosmetics, and environmental protection [START_REF] Valle | Cyclodextrins and their uses: a review[END_REF]. The cyclodextrins are cyclic oligosaccharides consisting of glucopyranose units (six α-cyclodextrin, seven β-cyclodextrin, eight γ-cyclodextrin, or more) linked by α-(1 → 4) bond, with a truncated conic shape [START_REF] Valle | Cyclodextrins and their uses: a review[END_REF]. The cyclodextrins have an inner cavity with a hydrophobic character, while the external structure has a hydrophilic character. During the inclusion, the guest molecules are held within the cavity of cyclodextrin, and this affects the physical and chemical properties of guests (solubility, stability…) [START_REF] Valle | Cyclodextrins and their uses: a review[END_REF][START_REF] Blum | Vitamins, 11. Niacin (nicotinic acid, nicotinamide)[END_REF].

Nicotinic acid (NA) (pyridine-3-carboxylic acid, vitamin B5) is a vitamin of plant origin food and plays an essential role in the metabolism of proteins, fats, and carbohydrates [START_REF] Blum | Vitamins, 11. Niacin (nicotinic acid, nicotinamide)[END_REF]. Nicotinic acid is broadly used as a lipid-modifying drug [START_REF] Carlson | Nicotinic acid: the broad-spectrum lipid drug. A 50th-anniversary review[END_REF]. The absence of this vitamin leads to pellagra disease. The NA is slightly soluble in water [START_REF] Haynes | CRC Hasndbook of Chemistry and Physics[END_REF] and exists in water in neutral and zwitterion form depending on pH and temperature [START_REF] Gonçalves | Solubility of nicotinic acid in the water, ethanol, acetone, diethyl ether, acetonitrile, and dimethyl sulfoxide[END_REF]. To enhance the solubility and protect this vitamin from external effects such as oxidation and structure modification, many works studied the inclusion of nicotinic acid with β-cyclodextrin (β-CD) [START_REF] Terekhova | Volume and heat capacity studies to evidence interactions between cyclodextrins and nicotinic acid in water[END_REF][START_REF] Terekhova | 1H NMR and calorimetric study on the binding ability of cyclodextrins to isomeric pyridine carboxylic acids in aqueous solution[END_REF][START_REF] Saha | Study to explore the mechanism to form inclusion complexes of β-cyclodextrin with vitamin molecules[END_REF]. In recent years, different researchers are interested [START_REF] Safia | Density functional theories study of the interactions between host β-Cyclodextrin and guest 8-Anilinonaphthalene-1-sulfonate: Molecular structure, HOMO, LUMO, NBO, QTAIM and NMR analyses[END_REF][START_REF] Gamboa-Carballo | Theoretical evaluation of the molecular inclusion process between chlordecone and cyclodextrins: A new method for mitigating the basis set superposition error in the case of an implicit solvation model[END_REF] in the detailed description of β-cyclodextrin host-guest systems and inclusion driving forces using quantum mechanical methods such as semi-empirical and density functional theory (DFT) [START_REF] Hadjer | Theoretical Investigation on Inclusion Complex of (s)-2-Isopropyl-1-(o-nitrophenyl) Sulfonyl) Aziridine with β-Cyclodextrin[END_REF]12]. Guendouzi et al. [START_REF] Guendouzi | Conformational Investigation of the Encapsulation of Nicotinic Acid Into β-Cyclodextrin[END_REF] investigate 1:1 encapsulation of NA with β-CD in the light of DFT using the Becke 97 method with D3 dispersion correction [START_REF] Grimme | Effect of the damping function in dispersion corrected density functional theory[END_REF] and inspect complex proprieties taking into account the ionic form of NA.

In the present work, we examine the 1:1 encapsulation of β-CD and the neutral form of nicotinic acid (Fig. 1) and describe the different driving forces using for the first time in the β-CD system the most recently global-hybrid meta-NGA of the Minnesota family the so-called MN-15 [START_REF] Haoyu | MN15: A Kohn-Sham global-hybrid exchange-correlation density functional with broad accuracy for multi-reference and single-reference systems and non-covalent interactions[END_REF] and hybrid meta-GGA the so-called PW6B95-D3 [START_REF] Zhao | Design of density functionals that are broadly accurate for thermochemistry, thermochemical kinetics, and nonbonded interactions[END_REF] (Perdew-Wang-91, 6-parameter functional exchange, and Becke-95 correlation with D3 Grimme's dispersion). The GGA exchange-correlation B3LYP [START_REF] Stephens | Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields[END_REF] (Becke, 3-parameter functional exchange, and Lee-Yang-Parr correlation) is also used as a benchmark with polarized Pople basis set 6-31G(d) in gas also in aqueous phase using the polarizable continuum model (PCM) [START_REF] Mennucci | Polarizable continuum model[END_REF].

The high level of theory was used after conformational research with different semi-empirical methods to define the global minimum structures of complexes. Moreover, the structural parameters were analyzed using the Akaike criterion [START_REF] Akpa | Small-sample corrected Akaike information criterion: an appropriate statistical tool for ranking of adsorption isotherm models[END_REF]. The frontier molecular orbital analyses are used to describe the charge flow and indicate the stability of the complexes. Thereafter, the dipole moment is used to confirm the orientation of the guest inside the β-cyclodextrin. In the population analysis, natural bond orbitals (NBO) are estimated. The quantum theory atom in molecule analysis (QTAIM) and non-covalent interaction analysis (NCI) were used to explain in detail the different interactions between the host and the guest.

The initial structures of β-cyclodextrin and nicotinic acid were downloaded from the PubChem database [START_REF] Kim | PubChem in 2021: new data content and improved web interfaces[END_REF], and the docking transition and rotation were achieved using the Hyperchem 7.51 software [START_REF] Hyperchem | 51 for windows 2002 Hypercube. Inc[END_REF]. The semi-empirical calculations of the conformational research were performed by MOPAC2016 [START_REF] Stewart | Stewart Computational Chemistry[END_REF], and the input files were generated by Avogadro [START_REF] Hanwell | Avogadro: an advanced semantic chemical editor, visualization, and analysis platform[END_REF], while the high level of theoretical calculations for the initial structure optimization, the final re-optimization, single point energies, and the NBO calculation were made by the Gaussian 16 software [24]. The results of visualizations were done by GaussView 5 [25]. Moreover, the QTAIM and NCI-RDG analysis was done by the Multiwfn [26] software and visualized by VMD [START_REF] Humphrey | VMD: visual molecular dynamics[END_REF]. Finally, Fig. 1 The β-cyclodextrin (host) and nicotinic acid (guest) structures Minitab2019 [START_REF]Minitab[END_REF] was used to plot the graphs and achieved the statistical calculations.

The identification of the most stable structure of the complex NA@β-CD (NA at β-CD) was described using manual conformational research [START_REF] Song | PM3 calculations on the complexation of α-cyclodextrin with the ground and excited quinone[END_REF] to get the potential energy surface. The glycosidic oxygen atoms of the β-CD were placed in the XY plane and their center was set up as a coordination system center; the wide rim was placed toward the Z-axis in a positive direction.

The nicotinic acid molecule enters through the wide rim of β-CD and vice versa, as two possible orientations were taken. The first orientation is from wide to narrow rim, namely orientation A, while the other one is from narrow to wide rim, namely orientation B (Fig. 2).

The bond between the C4 and C9 atoms of the guest was lined on the Z-axis, while the C4 atom is used as a center of the NA molecule. Then, the NA was moved in both orientations from -7 Å to 7 Å by 1 Å. In each step, the guest molecule was rotated around the Z-axis by 10° from 0° to 360°.

All the structures of NA@β-CD were fully optimized by three semi-empirical parametric methods: PM7 [START_REF] Stewart | Optimization of parameters for semiempirical methods VI: more modifications to the NDDO approximations and re-optimization of parameters[END_REF], PM6-D4H3 [START_REF] Rezac | Advanced corrections of hydrogen bonding and dispersion for semiempirical quantum mechanical methods[END_REF] (the PM6 method with D4 Grimme's dispersion and four hydrogen-bonding corrections), and PM3 [START_REF] Stewart | Optimization of parameters for semiempirical methods I[END_REF]. The complexation energy change implicated in the inclusion process produces two curves for two orientations for each method. The global minimum of both orientations given by the best method was re-optimized using the mentioned density functional theory (DFT) methods: MN-15, B3LYP, and PW6B95-D3 with 6-31G(d) in gas and aqua.

Results and discussion

The energy of complex

The formation energies of inclusion complex NA@β-CD which come from the low level of theory were calculated with the equation below [START_REF] Ohashi | Molecular mechanics studies on inclusion compounds of cyanine dye monomers and dimers in cyclodextrin cavities[END_REF]:

where E NA@β-CD is the optimized energy of complex, and E β-CD and E NA are the optimized energy of free host and guest respectively; the results are represented in Fig. 3.

As can be seen in Fig. 3, the PM7 gives a high level of energy in both orientations when compared with the other methods and it is worth mentioning that NA was located outside of β-CD in the global minima in both orientations, while the other methods give global minima, which NA embedded in the β-CD cavity, at 0 Å (250°) in case of PM6-D3H4 and -1 Å (250°) in case of PM3, in orientations A and B respectively (Table 1). The lower value energy is

(1) ΔE comp = E NA@β-CD -E NA + E β-CD

Fig. 2 A graphic representation of the inclusion process

Fig. 3 The potential energy surface for both orientations related to the most stable complex; PM3 gives the lower energy when compared with PM6-D3H4 with a difference in energy between the two minimums of 0.01 eV.

The lowest energy was given by PM3 and the two complexes of both orientations were re-optimized using MN-15\6-31G(d), B3LYP\6-31G(d), and PW6B95-D3\6-31G(d) level of theory in gas and aqueous phase.

The complexation energy of models A and B in the high level of theory can be calculated by Eq. ( 1).

The deformation energy describes the change of energies during the complexation of the separated host and guest, which were quantified by the above equation [START_REF] Barbiric | A molecular mechanics study of 1: 1 complex between azobenzene derivatives and β-cyclodextrin[END_REF][START_REF] Šponer | Computational approaches to the studies of the interactions of nucleic acid bases[END_REF]:

where E i (β-CD) and E i (NA) are respectively the single point energies of the isolated host and guest from the complex, and E β-CD and E NA are the optimized energies of the free host and guest respectively.

To identify the interaction energy during the complexation, the equation above was used [START_REF] Ohashi | Molecular mechanics studies on inclusion compounds of cyanine dye monomers and dimers in cyclodextrin cavities[END_REF]:

where E NA@β-CD is the optimized energy of the complex, and E i (β-CD) and E i (NA) are the single point energies of isolated host and guest from the complex respectively.

Table 2 displays the different values of energies for orientations A and B in gas and water. In the gas phase, (2)

ΔE def = (E i(β-CD) -E β-CD ) + (E i(NA) -E NA ) (3) ΔE int = E NA@β-CD -(E i(β-CD) + E i(NA) )
PW6B95-D3\6-31G(d) showed a lower value of complexation energy -33.8 kcal/mol and -27.32 kcal/mol in models A and B respectively; in the water phase, the complexation energy was -30.11 kcal/mol and -24.59 kcal/mol, in models A and B respectively.

Then, MN-15\6-31G(d) provides -30.79 kcal/ mol and -21.42 kcal/mol in gas, and -26.75 kcal/mol and -20.19 kcal/mol in aqua in orientations A and B respectively, whereas B3LYP\6-31G(d) reported significantly high complexation energy than the other methods. These results indicate that the most favorable orientation is orientation A, where the NA molecule is completely included from the carboxylic side toward the narrow rim of the host (Fig. 4).

When looking at the interaction energies, model B is higher in energy than A in MN-15\6-31G(d) about 9.51 kcal/ mol and 7.27 kcal/mol in gas and water respectively. However, B3LYP\6-31G(d) is higher about 8.47 kcal/mol and 4.95 kcal/mol in gas and water respectively. While PW6B95-D3\6-31G(d) showed opposite results, the orientation A is higher than B with a relatively small difference of 2.56 kcal/ mol and 3.3 kcal/mol in gas and water respectively.

The deformation energy was much higher in orientation A; this indicates the change of the structures and the gain of energy from the interaction during the complexation in both phases in MN-15\6-31G(d) and B3LYP\6-31G(d); however, in PW6B95-D3\6-31G(d), the orientation B was more deformed than A with the difference of ⁓ 9 kcal/mol in gas and water. These results favored orientation A in the case of MN-15\6-31G(d) and B3LYP\6-31G(d), while the orientation B is energetically favored in the case of PW6B95-D3\6-31G(d), and this shows that the results given by PW6B95-D3 are the opposite of the experimental data [START_REF] Terekhova | Volume and heat capacity studies to evidence interactions between cyclodextrins and nicotinic acid in water[END_REF][START_REF] Terekhova | 1H NMR and calorimetric study on the binding ability of cyclodextrins to isomeric pyridine carboxylic acids in aqueous solution[END_REF].

Geometry optimization

The comparison between the calculated and experimental bond lengths, bond angles, and the important dihedrals of the host and the guest were expected by linear regressing (Table 3).

The linear regressing exhibits a good correlation between experimental and calculated structural parameters of free β-CD and NA in the gas phase [START_REF] Boto | Revealing strong interactions with the reduced density gradient: a benchmark for covalent, ionic and charge-shift bonds[END_REF][START_REF] Lane | Are bond critical points critical for hydrogen bonding?[END_REF]. The coefficient of determination R 2 of bond lengths model is over 98% and 97% for NA and primary hydroxyl bonds of β-CD; it is worth mentioning that the primary hydroxyls were selected because they interfere in the complexation process. On the other hand, R 2 of the angle model of the nicotinic acid molecule is between 94.9 and 95.5%, while in the case of β-cyclodextrin, R 2 of angles and dihedrals were ranged between 8.24 and 9.8%.

The Akaike information criterion correction (AICc) has given more insight into the comparison of the given results of each method. AICc was used because of its strong ability to operate with small available data taken from the experiment in the case of NA as well as the limited structural parameter selected from the β-CD molecule [START_REF] Khan | Density functional theory study of emerging pollutants removal from water by covalent triazine based framework[END_REF]. The best fit model is the lowest AICc value (Table 3).

The ranking of structure models is as follows: MN-15 gives a good correlation in the case of bond length, secondly, PW6B95-D3, and at least B3LYP in NA and β-CD. While in the case of bond angles, B3LYP had the smallest AICc and should be the first for NA structure, secondly, MN-15, and the last one is PW6B95-D3, unlike the host AICc of MN-15 and PW6B95-D3 respectively were smaller than the B3LYP method and this ranked MN-15 method firstly.

The calculated dihedral angles of the β-cyclodextrin showed high AICc and small R 2 approximately equal to 8%; this revealed a bad correlation with the experiment at the high level of theory.

The regression of structural parameters between free and complex molecules in water (Table 4) indicated that the β-CD primary hydroxyl bond lengths and bond angles in orientation B changed more than their NA in MN-15\6-31G(d); in contrast, the NA structural parameter at B3LYP\6-31G(d) and PW6B95-D3\6-31G(d) was changed in orientation A than in B. However, the host selected dihedral angles changed in A more than B. In addition, the selected dihedrals of NA make a significant change in both models. These reveal that the inclusion occurred and the β-CD primary hydroxyls play a significant role in the inclusion process.

The synthesis of complex NA@β-CD originally occurs in water; as a result, the calculation in the aqueous phase was considered in the rest of this work.

The frontier molecular orbital analysis

To further understand the inclusion process and to determine the stability of the complex in the light of density function theory, the analysis of frontier orbitals has been used. The highest occupied molecular orbital HOMO represents the ability of the molecule to give an electron, while the lowest unoccupied molecular orbital LUMO represents the ability of the molecule to accept an electron. The energies of HOMO and LUMO are related to ionization potential (I) and electronic affinity (A) respectively, according to Koopman's theory [START_REF] Koopmans | Über die Zuordnung von Wellenfunktionen und Eigenwerten zu den einzelnen Elektronen eines Atoms[END_REF] which:

where E HOMO and E LUMO are the HOMO and LUMO energy respectively. The HOMO is located in β-cyclodextrin while the LUMO is located in nicotinic acid for both models A and B. Figure 5 shows the frontier orbitals of the complex in the high level of theory.

The energy gap (E gap ) between HOMO and LUMO is an important factor that helps to determine the stability of the system. The system with the higher gap is the most stable [START_REF] Pearson | Ionization potentials and electron affinities in aqueous solution[END_REF].

The global reactivity descriptors as electronic chemical potential (μ) [38, 39, 40, 41, ], global hardness (η) [START_REF] Pearson | Absolute electronegativity and hardness: applications to organic chemistry[END_REF], global electrophilicity (ω) [START_REF] Parr | Electrophilicity index[END_REF], and electrophilicity-based charge transfer (ECT) [START_REF] Padmanabhan | Electrophilicity-based charge transfer descriptor[END_REF] were obtained at the base of molecular frontier orbitals:

where ΔN max is the maximum amount of electronic charge:

(4) E HOMO ≈ -I (5) E LUMO ≈ -A (6) = - I + A 2 ≈ E LUMO + E HOMO 2 (7) = I -A 2 ≈ E LUMO -E HOMO 2 (8) = 2 2 (9) ECT = ΔN max (host) -ΔN max (guest) (10) ΔN max = -
Table 5 shows frontier molecular orbital energies, the gap between HOMO and LUMO energies, and the calculated reactivity descriptors of the free β-cyclodextrin, free nicotinic acid, and the complex of orientations A and B.

In the encapsulation process, the electrons flow from the molecule with the high electronic chemical potential to the molecule with the low electronic chemical potential, till the equilibrium state. In our case, the charge is transferred from the host to the guest molecule. The negative value of the chemical potential of the complexes indicates the stability of the system.

The hardness measures the resistance to change in electron distribution in the molecule, the complex with the larger gap energy has a high hardness, and therefore the complex is more stable. The complex in orientation A is more stable than B because orientation A is harder.

The electrophilicity index is a measure of the power to attract electrons, the guest molecule is more electrophilic than the host, and also complex A was more electrophilic than B.

To confirm the charge transfer direction, the electrophilicity-based charge transfer, ECT, has been used. The value of ECT was negative, which means the charge flows from β-CD to NA.

We conclude from the previous results that orientation A at all DFT methods used in this work-including PW6B95-D3\6-31G(d)-is the favored complex.

The Cartesian and spherical coordination of the dipole moments is listed in Table 6; the dipole moment describes the asymmetric charge distribution of the molecules and describes the orientation of the complex in the light of dipole-dipole interaction.

The β-CD dipole moments in orientation A and orientation B are approximately within 7 Deby to 11 Deby. The value of the dipole moment is lower in orientation A according to the three DFT methods; in addition, the dipole moment of free molecules changed during the inclusion process due to the polarization [START_REF] Hamai | Inclusion effects of cyclomaltohexa-and heptose (α-and β-cyclodextrins) on the acidities of several phenol derivatives[END_REF]. The antiparallel direction of the dipole moments of host and guest played an essential role to reveal the guest orientation and the stability of the complex [START_REF] Kitagawa | A Molecular Orbital Study of Cyclodextrin Inclusion Complexes. I. The Calculation of the Dipole Moments of α-Cyclodextrin Aromatic Guest Complexes[END_REF].

Table 4 The regression coefficient of structural change of nicotinic acid and β-cyclodextrin primary hydroxyls in orientations A and B in the aqueous phase a Bond length; b bond angles; c dihedral angles To describe the contribution of the dipole moments in the encapsulation process, the Z-axis dipole moment (μ Z ) is used because the major contribution of the vectors was along the Z-axis than x and y (Table 6), whether of isolated molecules or the complexes. In addition, the azimuthal angle φ can help to indicate how parallel the Z-axis dipolar moment is. A and B respectively. On the other hand, Z-axis dipolar moments of NA at MN-15\6-31G(d) and PW6B95-D3\6-31G(d) in the orientation A are located away from the Z-axis in the negative direction by 24.3° and 25.6° respectively, and in the orientation B, the μ Z is located away by 69.6° in the negative direction but in the case of the PW6B95-D3\6-31G(d) and B3LYP\6-31G(d), the angles φ are 87° and 83.8° respectively. These results emphasize clearly that the dipole moments in orientation A were more antiparallel than those of orientation B. Moreover, the dipolar moments in orientation A were smaller than those in B in the selected DFT methods. This describes the stability of the complex due to the dipole-dipole interaction between host and guest, as a driving force of the inclusion

Orientation A Orientation B R 2a (%) R 2b (%) R 2c (%) R 2a (%) R 2b (%) R 2c (%) β-CD MN-15\6-31G(

Natural population analysis

The charge distribution within the nicotinic acid before and after the complexation was achieved using the natural population analysis (NPA) [START_REF] Reed | Natural population analysis[END_REF] and is listed in Table 7. This method has an important role in quantum chemistry applications which effectively detects reactive sites and other molecular interactions including the charge transfer by treating the natural bond orbitals [START_REF] Gupta | Principles and applications of quantum chemistry[END_REF], unlike the Mulliken population which treats all orbitals and is affected by the basis set [START_REF] Wiberg | Comparison of atomic charges derived via different procedures[END_REF].

What stands out in Table 7 is the negative charge of the nicotinic acid molecule is located mainly in the oxygen of the carboxylic function group: O149 of carbonyl with -0.71 and the O148 attached to the acidic hydrogen with -0.63, also in the pyridine ring which the negative charge ranged from -0.14 to -0.48, are located respectively in N150, C154, C151, and C152. On the other hand, the positive charge ranging from 0.04 to 0.83 is located in the rest of the atoms: C156, all hydrogens (especially in H161 which has the maximum magnitude of 0.52), C153, and C155 respectively.

During the inclusion, the charge distribution changed in both models; the change was estimated by the difference between the charge of the molecule in the complex and the free molecule given by the mentioned methods; Fig. 6 shows these results.

The change in orientation A is relatively superposed in both methods where the O149, O148, and C152 exhibit increasing in the negative charge which acted as an acceptor; in addition, the most positive change was exhibited in the carbonyl carbon C156 and increased in acidic hydrogen H161, unlike the other hydrogens. However, the orientation B gives a noticeable contrast in charge changing between the methods, such as the difference between the increasing or the decreasing of charge in some atoms, in addition to the opposite behavior of some atoms, like O149 which acts as an acceptor in PW6B95-D3\6-31G(d) and as a donor in MN-15\6-31G(d) and vice versa in the case of H161. In general, the charge transfer takes place from the β-CD to the NA in orientation A where it reveals a relatively high negative value unlike orientation B where the charge transfer takes place from the NA to the β-CD and gives a relatively small positive value, and this favored the orientation A more than B, and therefore this supports the results mentioned above.

Natural bond orbital analysis

Natural bond orbital (NBO) analysis emphasizes the role of intermolecular orbital interactions, especially charge transfer. The delocalization of orbitals can be explained in terms of stabilization energy and the interaction intensity between the occupied (donor) and unoccupied (acceptor) orbitals are directly related to its value [START_REF] Reed | Intermolecular interactions from a natural bond orbital, donor-acceptor viewpoint[END_REF]. The energy is estimated by the second-order perturbation theory E (2) : where q i is the donor orbital occupancy, ε i and ε j are diagonal elements (orbital energies) and F (i,j) is the offdiagonal NBO Fock matrix element. The energies are listed in Table 8.

In orientation A, the delocalized orbitals are the same in MN-15\6-31G(d) and PW6B95-D3\6-31G(d), while orientation B gives a different delocalization for each method. The main contribution of β-CD is acting as a donor, and the primary hydroxyl oxygen lone pairs are involved. The stabilization energy is much higher when the charge transfer occurs from the lone pair of O11 and the anti-bonding orbital O148-H161, the lone pair of O33 also involved.

When the guest acts as a donor and the host as an acceptor, the transfer of charge arises especially from O149 and N150 lone pairs to σ * orbitals of the host according to the values of the stabilization energy. It is worth mentioning that the stabilization energies given by the two methods in the case of orientation A are too close, and some of these values can confirm the existence of hydrogen bonding.

The orientation B gives different orbital delocalization depending on the chosen method. MN-15\6-31G(d) shows low stabilization energy compared with orientation A from the same method, while PW6B95-D3\6-31G(d) gives higher energy and mostly arises when the β-CD acts as a donor, especially from the lone pair of secondary hydroxyl oxygen O71 to the anti-bonding orbital of O148-H161 in NA; on the other hand, when NA acts as a donor, the charge transfer occurs from the lone pair of N150 to the anti-bonding orbital of O22-H97. The results confirm that the charge transfer drives the stability of the host-guest system. ( 11)

E (2) = q i F (i,j) j -i

Quantum atom in molecule analysis

The main driving force in the inclusion system is the noncovalent interactions between the host and the guest; these interactions lead to a stable complex. Bader's quantum theory atom in molecule (QTAIM) [START_REF] Bader | Atoms in molecules-a quantum theory[END_REF] gives a crucial description of the interactions in the molecular system, based on the topology analysis of electronic density and its reduced gradient at the critical points. We note that the critical point is the point where the gradient of electronic density vanishes [START_REF] Nkungli | Theoretical analysis of the binding of iron (III) protoporphyrin IX to 4-methoxyacetophenone thiosemicarbazone via DFT-D3, MEP, QTAIM, NCI, ELF, and LOL studies[END_REF]. We are interested in this study in the bond critical points (BCP). Whereas the Hessian matrix of BCP is characterized by three non-zero eigenvalues, two are negative and one is positive. The electronic density ρ b and its corresponding Laplacian ∇ 2 ρ b provide information about the nature and strength of bonding between host and guest. The covalent bonding is characterized by a large ρ b value and ∇ 2 ρ b < 0, whereas the small ρ values and ∇ 2 ρ b > 0 character- ize the closed-shell interactions [START_REF] Bader | The characterization of atomic interactions[END_REF]. Furthermore, Cremer and Karka [START_REF] Cremer | Theoretical determination of molecular structure and conformation. 15. Three-membered rings: bent bonds, ring strain, and surface delocalization[END_REF] found that the total electronic energy density H b (RC) at the bond critical point helps to characterize the type of bonding. It is worth adding in the context that the total electronic energy density H b (RC) is the sum of kinetic energy densities G b (RC) (always positive) and potential energy densities V b (RC) (always negative) at BCPs which [START_REF] Cremer | Theoretical determination of molecular structure and conformation. 15. Three-membered rings: bent bonds, ring strain, and surface delocalization[END_REF]. The total electronic energy density H b (RC) has a negative value when the concentration of electronic charge is large in the inter-nuclear region (it means G b (RC) < -V b (RC)), and this indicates a covalent bond or strong hydrogen bond [START_REF] Rozas | Behavior of ylides containing N, O, and C atoms as hydrogen bond acceptors[END_REF]; on the other hand, the total energy has a positive value, when the charge depletes away from the inter-atomic region to the nuclei (it means G b (RC) > -V b (RC)), and this indicates closed-shell interactions such as medium and weak hydrogen bond including van der Waals interactions. Nakanishi and coworkers developed a criterion [START_REF] Nakanishi | Atoms-in-molecules dual parameter analysis of weak to strong interactions: behaviours of electronic energy densities versus laplacian of electron densities at bond critical points[END_REF][START_REF] Nakanishi | Polar Coordinate Representation of H b (r c ) versus (ℏ 2 /8 m)▽ 2 ρ b (r c ) at BCP in AIM Analysis: Classification and Evaluation of Weak to Strong Interactions[END_REF] which enables to analyze the type of interactions based on the electronic density ρ b , Laplacian of the electronic density ∇ 2 ρ b , and the total electronic energy density H b (RC), where these parameters should lie in the range of 0.01 au < ρ b < 0.04 au, 0.06 au < ∇ 2 ρ b < 0.12 au, and -0.01 au < H b (RC) < 0.003 au; in the case of hydrogen bonding, as well in the case of van der Waals interactions, the range should be 0.00 au < ρ b < 0.01 au, 0.00 au < ∇ 2 ρ b < 0.05 au, and 0 au < H b (RC) < 0.003 au. Espinosa et al. [START_REF] Espinosa | Hydrogen bond strengths revealed by topological analyses of experimentally observed electron densities[END_REF] formulated a proportionality between the interaction energy and the potential energy at the bond critical point where E int = 0.5V b (RC), which is helpful to classify the non-covalent interaction, especially the hydrogen bonds. In terms of interaction energy, the hydrogen bonds are estimated to be strong H-bond (15-45 kcal/mol), medium H-bond (4-15 kcal/mol), and weak H-bond (1-4 kcal/mol) [START_REF] Gilli | The nature of the hydrogen bond: outline of a comprehensive hydrogen bond theory[END_REF]. Regarding the van der Walls interactions, the energy varies in the range of 0.1-4 kcal/mol [START_REF] Atkins | Atkins' physical chemistry[END_REF].

H b (RC) = G b (RC) + V b (RC)
Tables S14, S15, S16, and S17 (supplementary data) summarize the topological parameters of both orientations A and B of the two methods MN-15\6-31G(d) and PW6B95-D3\6-31G(d) in the aqueous phase that has been selected, named ρ b ,∇ 2 ρ b , H b (RC), G b (RC), V b (RC), and E int . The interact atoms mentioned in the table and the type of interaction thereof. The bond critical points and the bond paths are provided in Fig. 7.

According to the data listed in the tables mentioned above, all the interactions between β-CD and NA in all cases are noncovalent in nature, since the corresponding values of ρ b are small and ∇ As can be seen in the tables, the interactions present between β-CD and NA are O

•••H, N•••H, C•••H, = O•••H, O•••O, N•••O,

Non-covalent interaction-reduced density gradient analysis

The non-covalent interaction method (NCI) is the most frequent tool used for visualizing the non-covalent interactions in real space [START_REF] Johnson | Revealing noncovalent interactions[END_REF] such as the hydrogen bonds, van der Waals interactions, and steric effects in light of the reduced density gradient (RDG), given by this equation: This quantity is dimensionless and has been used as a bonding descriptor that describes the deviation of homogeneous density distribution in the molecular system; the RDG has a very low value at bond critical points [START_REF] Boto | Revealing strong interactions with the reduced density gradient: a benchmark for covalent, ionic and charge-shift bonds[END_REF] that could be associated with the presence of covalent bonds and noncovalent interactions; conversely, RDG in regions around the nuclei or far from the molecule takes a large positive value; the RDG isosurface enables to reveal the nature and strength of the non-covalent interactions using the sign of second eigenvalues of electron density Hessian matrix λ 2 multiplied by the electronic density ρ (Sign(λ 2 )ρ), that for the attractive interactions Sign(λ 2 )ρ < 0, and the repulsive interactions Sign(λ 2 )ρ > 0; these are because of the negative value of λ 2 associated with the bond critical points and on the other hand the positive value of λ 2 associated with the ring critical points respectively [START_REF] Lane | Are bond critical points critical for hydrogen bonding?[END_REF].

The 3D visualization of NCI interaction is mapped using Sign(λ 2 )ρ to show the RDG isosurface in real space, with a color-coding scheme as follows: blue-green-red which represents respectively H-bonds or halogen bonds, van der Waals interactions, and steric repulsions. Moreover, the RDG vs Sign(λ 2 )ρ can be represented in 2D scatter diagram for all the interactions of the molecular system [START_REF] Khan | Density functional theory study of emerging pollutants removal from water by covalent triazine based framework[END_REF].

As can be seen in orientation A in MN-15\6-31G(d) and PW6B95-D3\6-31G(d) level of theory represented in Fig. 8, which highlighted the different NCI between the host and the guest, the O11 of the host primary hydroxyl groups and O149 oxygen of the guest carboxylic functional group correspond to strong and medium hydrogen bonds that appear in blue and blue-green disk shape by respect. The weaker interaction van der Waals is observed as a flat green isosurface mainly around the upper of the pyridine ring related to H-H bonding and π-H-bonding. In orientation B, the N150 Fig. 8 The down view of 3D-RDG of the complex A and complex B in water given by MN-15\6-31G(d) (up) and PW6B95-D3\6-31G(d) corresponds to a strong hydrogen bond which appears in blue disk shape while the weak interactions as a green scattered surface are mainly located around the carboxylic functional group also distributed around the hydrogens of the ring in the case of MN-15\6-31G(d). On the other hand, PW6B95-D3\6-31G(d) showed strong to medium hydrogen bonds appearing in blue to blue-green disk-shaped concerns to O11, N150, and O149 respectively; also a green scattered surface around the hydrogens of the NA molecule corresponds to the van der Waals interaction. In orientations A and B in both methods, half of the disk shape in green indicates a weak attractive interaction between the hydrogen atoms of the pyridine ring and the carboxylic oxygens, and the other half in red expects a steric repulsion between the carbons and oxygens of the same molecule. In addition, a tiny red oval-shaped surface is remarked inside the nicotinic acid ring indicating steric repulsion.

In the 2D plot in Fig. 9, the stronger interactions are seen as narrow spikes to the left with higher density values ρ > 0.01 au which are represented in blue; in the middle of Fig. 9 The down view of 2D-RDG of the complex A and complex B in water given by MN-15\6-31G(d) and PW6B95-D3\6-31G (d) the plot where the ρ < 0.01 au, the van der Walls are illustrated in green, while the steric repulsion represented in the right side in red. A strong to medium H-bonds showing as patched blue dots correspond to O11•••H161 within -0.03 to -0.05 au, and blue spikes correspond to O148•••H147 at -0.02 to -0.03 au, in the orientation A given by the two methods. While in the case of MN-15\6-31G(d), a broad spick appears at -0.03 au corresponding to N150•••H97 as medium H-bond, but in PW6B95-D3\6-31G(d), blue spikes patched within -0.045 au < Sign(λ 2 )ρ < -0. The spikes in green at Sign(λ 2 )ρ > -0.01 au correspond mainly to the van der Waals interactions between the ring of NA and the host cavity in the case of orientation A; on the other hand, in the case of orientation B, this region corresponds to weak interactions acts mainly at the hydrogen atoms of the guest or the host. Also, it can be noticed as well a spike at 0.02 au corresponds to the ring steric repulsion in the red region in both orientations and methods.

These results demonstrate the presence of non-covalent interactions between the β-CD and NA as a driving force and elucidate the important role of strong interactions in controlling the guest orientation and the role of the weak interactions in the stability of the complexes.

Conclusion

A computational investigation on the inclusion of 1:1 stoichiometry of β-cyclodextrin and nicotinic acid in neutral form has occurred using different quantum chemistry methods of calculation and analysis. The conformational research was carried out using different parametric semi-empirical methods such as PM7, PM6-D3H4, and PM3. The most stable conformations were re-identified by DFT Minnesota family method MN-15, B3LYP, and three dispersion-correction PW6B95 using polarized Pople basis set 6-31G(d). Orientation A gives the favored complex, in which the guest was encapsulated and the carboxylic group of nicotinic acid was oriented toward the narrow rim of β-cyclodextrin; the orientation B shows fewer interaction energies in the case of PW6B95-D3\6-31G(d). However, the structural parameters described the structure deformation of the host and the guest during the inclusion and confirmed the important role of the primary hydroxyl rim of β-cyclodextrin in the inclusion process. Also, orientation A gives the harder complex which has the higher gap values. In addition, the dipole-dipole interaction favored orientation A and helped to identify the guest orientation. The polarization of the host and the guest are changed during the inclusion and the charge transfer was carried out from the host to the guest in all mentioned DFT methods. Finally, the hydrogen bonds are observed in the complexes at MN-15\6-31G(d) and PW6B95-D3\6-31G(d), and the van der Waals interactions exist.

Fig. 4

 4 Fig.[START_REF] Haynes | CRC Hasndbook of Chemistry and Physics[END_REF] The complex of orientation A (left) and B (right) in MN-15\6-31G(d) level of theory

  The β-CD azimuthal angles at MN-15\6-31G(d) and PW6B95-D3\6-31G(d) in orientation A are 43.9° and 41.5° respectively, while in the orientation B, the azimuthal angles are 29° and 36.2° respectively. But at B3LYP\6-31G(d) level of theory, μ Z is far from the Z-axis in the opposite direction by 22.4° and 44.2° in orientations

Fig. 5

 5 Fig. 5 A graphic representation of frontier molecular orbitals in the level of theory in water

Fig. 6

 6 Fig. 6 The change of charges after the inclusion of nicotinic acid in MN-15\6-31G(d) and PW6B95-D3\6-31G(d) in orientations A and B in the aqueous phase

  2 ρ b are all positive, in orientation A within 0.0039 au ≤ ρ b ≤ 0.045 au and 0.0026 au ≤ ρ b ≤ 0.045 au in MN-15\6-31G(d) and PW6B95-D3\6-31G(d) respectively, and, in the same way, 0.015 au ≤ ∇ 2 ρ b ≤ 0.14 au and 0.011 au ≤ ∇ 2 ρ b ≤ 0.14 au for MN-15\6-31G(d) and PW6B95-D3\6-31G(d) respectively. On the other hand, the values in orientation B fall within 0.0033 au ≤ ρ b ≤ 0.034 au and 0.002 ≤ ρ b ≤ 0.043 au according to MN-15\6-31G(d) and PW6B95-D3\6-31G(d) respectively. Also, the Laplacian ranged in 0.013 au ≤ ∇ 2 ρ b ≤ 0.094 au and 0.0087 au ≤ ∇ 2 ρ b ≤ 0.13 au for MN-15\6-31G(d) and PW6B95-D3\6-31G(d) respectively. Furthermore, the energy at the critical point H b (RC) shows the existence of two kinds of bonding:H-bonds and van der Waals interaction, ranging in orientation A as -0.0022 to 0.0017 au and -0.0023 to 0.0016 au, using MN-15\6-31G(d) and PW6B95-D3\6-31G(d) respectively, while in orientation B, the energy at BCP lies within -0.0015 to 0.0017 au and -0.0025 to 0.0019 au in MN-15\6-31G(d) and PW6B95-D3\6-31G(d) respectively.

  and H•••H; there are no differences in the interacted atoms between the structures obtained from the two methods in the complex A, except H(157)•••H(144) which is given by MN-15\6-31G(d). However, in complex B, the structures obtained by the DFT methods give several similar interacted atoms with different interaction types, and this depends on the deformation of the β-CD and NA structures. The β-cyclodextrin clasped the nicotinic acid molecule with the upper and the lowest hydroxyls in PW6B95-D3\6-31G(d) and this led to deforming the structure of NA, and this gives different interact atoms: N150•••H104,

Fig. 7

 7 Fig.7The bond critical points and the bond paths of the complex A and complex B in water given by MN-15\6-31G(d) and PW6B95-D3\6-31G(d) (down view)

  02 au correspond from strong to medium H-bonds: O71•••H161, N150•••H97, and O149•••H140 by respect. However, weak H-bonds have been observed within -0.02 au < Sign(λ 2 )ρ < -0.01 au in green, in orientation A at the two methods which correspond to O149•••H144 and N150•••H114. While in orientation B, the weak H-bonds can be noticed at the same region in the plot corresponding to O11•••H158, O149•••H134, and O148•••H84 at MN-15\6-31G(d) level of theory, and O149•••H141 at PW6B95-D3\6-31G(d).

Table 1 The

 1

	most stable conformations for each semi-	Orientation	PM7		PM6-D3H4		PM3
	empirical method		A	B	A	B	A	B
		Position (Å)	-6	-7	0	2	0	-1
		Angle (°)	40	230	250	20	170	250
		Energy (eV)	-0.29	-0.21	-0.70	-0.45	-0.69	-0.71

Table 3 The

 3

	Akaike information criterion correction		β-CD			NA
	and the regression coefficient of the nicotinic acid and β-cyclodextrin primary hydroxyls in the high level of theory in the gas phase	Bond length AICc R MN-15\6-31G(d) -74.8 97.9	Bond angle 60.9 8.2	Dihedral angles 179.6 8	Bond length -50.1 98.6	Bond angle 45.7 94.9
		B3LYP\6-31G(d)	-72.8 97.8	69.2 9.8	179.4 8.3	-49.3 98.5	45	95.5
		PW6B95-D3\6-31G(d) -73.4 97.8	64.4 9.4	179.5 8.2	-49.6 98.5	46	95.1

2 (%) AICc R 2 (%) AICc R 2 (%) AICc R 2 (%) AICc R 2 (%)

Table 5

 5 The global reactivity indexes of β-cyclodextrin and nicotinic acid in MN-15\6-31G(d), B3LYP\6-31G(d), and PW6B95-D3\6-31G(d) of free molecules and orientations A and B in the aqueous phase

		free β-CD		free NA		A			B
	MN-15\6-31G(d)/B3LYP\6-31G(d)/PW6B95-D3\6-31G(d)
	E HOMO (eV)	-0.295/ -0.257/ -0.270	-0.307/ -0.268/ -0.284	-0.296/ -0.255/ -0.271	-0.296/ -0.252/ -0.267
	E LUMO (eV)	0.084/0.052/0.062	-0.033/ -0.062/ -0.052	-0.038/ -0.063/ -0.055	-0.041/ -0.070/ -0.059
	E gap (eV)	0.380/0.308/0.332	0.275/0.207/0.232	0.258/0.191/0.216		0.255/0.182/0.208
	μ (eV)	-2.867/ -2.791/ -2.836	-4.623/ -4.489/ -4.562	-4.549/ -4.330/ -4.433	-4.582/ -4.378/ -4.440
	η (eV)	5.166/4.196/4.516	3.738/2.815/3.159	3.507/2.603/2.933		3.469/2.473/2.832
	ω (eV)	0.796/0.928/0.890	2.859/3.579/3.294	2.951/3.601/3.350		3.026/3.875/3.481
	ΔN max	0.555/0.665/0.628	1.237/1.595/1.444	/			/
	ECT	-0.682/ -0.930/ -0.816
	Table 6 The dipole moments in the Cartesian and spherical			μ X (Deby)	μ Y (Deby)	μ Z (Deby)	μ (Deby)	θ (°)	Φ (°)
	coordination in the aqueous phase of free β-cyclodextrin and nicotinic acid, and in orientations A and B, in addition to the complexes in A	MN-15\6-31G(d) Free β-CD Free NA β-CD in orientation A	-3.96 -0.10 -2.82	-1.82 0.66 4.09	1.96 0.00 5.17	4.78 0.66 7.17	-155.28 98.71 124.55	65.84 90.02 43.85
	and B		NA in orientation A	0.03	0.44	-0.96	1.05	85.77	155.35
			β-CD in orientation B	3.85	3.12	8.95	10.23	39.07	28.97
			NA in orientation B	-0.11	0.48	-0.18	0.52	102.62	110.36
			Complex A		-3.72	2.10	5.79	7.20	150.59	36.38
			Complex B		4.40	1.35	9.02	10.13	17.06	27.07
			B3LYP\6-31G(d)
			Free β-CD		1.85	-1.43	6.06	6.49	-37.78	21.10
			Free NA		-0.19	0.72	-0.00	0.75	104.79	90.02
			β-CD in orientation A	-3.16	-1.09	-8.10	8.76	-161.06	157.58
			NA in orientation A	0.05	0.63	1.33	1.47	84.74	25.31
			β-CD in orientation B	-1.37	-7.53	-7.87	10.98	-100.31	135.80
			NA in orientation B	-0.17	0.69	0.08	0.71	103.56	83.83
			Complex A		2.27	2.89	-9.15	9.86	51.84	158.12
			Complex B		3.63	6.05	-8.14	10.77	120.94	139.07
			PW6B95-D3\6-31G(d)
			Free β-CD		-3.66	-1.84	2.11	4.61	-153.30	62.73
			Free NA		-0.16	0.70	0.00	0.72	103.02	90.02
			β-CD in orientation A	-3.42	3.27	5.38	7.17	136.29	41.35
			NA in orientation A	0.03	0.46	-0.95	1.06	86.35	154.43
			β-CD in orientation B	2.73	5.13	7.93	9.83	61.97	36.22
			NA in orientation B	0.00	0.35	0.02	0.35	90.02	86.96
			Complex A		-3.94	1.17	6.01	7.28	163.45	34.36
			Complex B		-3.78	-4.12	8.04	9.79	-132.55	34.81

Table 7

 7 The natural charge population values of each atom of nico-

	tinic acid in MN-15\6-31G(d) and PW6B95-D3\6-31G(d) in free and
	orientations A and B in the aqueous phase
	Atoms	Free	A	B
	MN-15\6-31G(d)/PW6B95-D3\6-31G(d)
	O148	-0.711/ -0.712	-0.727/ -0.723	-0.721/ -0.724
	O149	-0.626/ -0.629	-0.677/ -0.68	-0.615/ -0.662
	N150	-0.478/ -0.479	-0.48/ -0.479	-0.519/ -0.516
	C151	-0.227/ -0.223	-0.214/ -0.211	-0.219/ -0.210
	C152	-0.147/ -0.152	-0.160/ -0.164	-0.142/ -0.147
	C153	0.065/0.061	0.068/0.06	0.072/0.081
	C154	-0.275/ -0.276	-0.271/ -0.272	-0.274/ -0.280
	C155	0.043/0.039	0.042/0.039	0.062/0.057
	C156	0.832/0.833	0.858/0.856	0.836/0.848
	H157	0.265/0.268	0.258/0.262	0.256/0.262
	H158	0.242/0.245	0.238/0.241	0.249/0.248
	H159	0.26/0.263	0.258/0.261	0.26/0.263
	H160	0.235/0.238	0.237/0.238	0.25/0.251
	H161	0.523/0.523	0.534/0.535	0.522/0.534
	Total	0/0	-0.035/ -0.033	0.017/0.004

Table 8

 8 Calculated NBO values for orientations A and B in MN-15\6-31G(d) and PW6B95-D3\6-31G(d) in the aqueous phase

	E (2) (kcal mol -1 )

Author contribution

The semi-empirical calculation, QTAIM and NCI-RDG studies, data collection, and analysis were performed by Houdhaifa R. Belhouchet. The final optimization using DFT and NBO analysis calculation was performed by Tahar Abbaz and Didier Villemin. The first draft of the manuscript was written by Houdhaifa R. Belhouchet and Tahar Abbaz; Amel Bendjedou commented on previous versions of the manuscript. Houdhaifa R. Belhouchet, Tahar Abbaz, Amel Bendjedou, and Abdelkrim Gouasmia read and approved the final manuscript.