Computational investigation of dimethoate and β-cyclodextrin inclusion complex: molecular structures, intermolecular interactions, and electronic analysis

The proposed study concerns the inclusion complexation of dimethoate (DMT) in the β-cyclodextrin (β-CD) molecule cage using a 1:1 stoichiometry. The interactions between DMT and β-CD were evaluated using PM7 and DFT in water and gas, using the CAMB3LYP functional 6-31G(d,p) basis set. All approaches agree with the optimal 3D structure, which includes full DMT inclusion in the β-CD cavity. Complexation, LUMO, and HOMO energies were computed. The natural bond orbital (NBO) and UV–visible calculations were determined and discussed. Additionally, the non-covalent intermolecular interactions between dimethoate and β-cyclodextrin are investigated through RDG, NCI, and IGM. The main forces stabilizing the examined inclusion complex are H-bond and van Der Waals interactions. Furthermore, the energy decomposition analysis (EDA) was applied highlighting the H-bonding by quantifying its strength. The TD-DFT method provided the electronic UV–vis spectra showing significant electronic transitions between the host and the guest by observing the molecular orbitals localizations.


Introduction
Recently, the scientific community certified that phytosanitary agents and their metabolites harm the environment. The use of pesticides in agriculture presents an ecotoxicological risk [1,2]. They found that pesticides regularly exceed ecological quality thresholds [3][4][5] and constitute a hazard for aquatic species and sediment. Besides, they threaten human health through water, fish, and marine species consumption [6][7][8][9][10].
Organophosphate pesticides are the most widely used insecticides, causing environmental degradation. In agriculture, the ecological problems of pesticides are created by high doses and their frequent applications [11].
Indeed, dimethoate is an organophosphate pesticide used to combat the losses due to pests like fungi and insects [12]. Inorganic insecticide dimethoate (Fig. 1a) potentially harms the nervous system. Both an insecticide and an acaricide serve their purpose well. Insects and human neurological systems depend on cholinesterase consisting of an enzyme these chemicals aim to prevent from doing their job.
Pesticide pollution of the aquatic ecosystem can result in acute and chronic intoxication of fish and other species.
Therefore, dimethoate is complexed with supramolecular structures involving inclusion complexes with cyclodextrins to remediate the pesticide pollution problem. Cyclodextrins are commonly used for water treatment. They retain a lot of pesticides, insecticides, metals, or toxic organic compounds such as phenols, thus purifying running water [16][17][18].
Cyclodextrin is a molecule cage of natural origin that creates an opportunity to encapsulate different compounds. Microcapsule technology's carrier matrix, cyclodextrin, and its derivatives are widely employed in many culinary and medicinal goods, prompting much scientific investigation [17].
Because of the particular molecular structure of cyclodextrin, which consists of an exterior hydrophilic surface and an interior hydrophobic cavity, they can form inclusion complexes with various guest molecules that have the appropriate polarity and dimensions [20][21][22][23][24][25].
As known, the β-CD (Fig. 1b) is the most commonly used for the formation of inclusion complexes with a variety of products through the interactions of van der Waals types (vdW) [26]. Cyclodextrin has an internal cage and a truncated cone with a diameter of 6-6.4 Å and a depth of 8 Å. It is important to develop inclusion complexes of organic compounds with cyclodextrins for technological and medicinal aims [27][28][29][30][31][32][33].
Some hydrophobic organic pesticides have limitations for extensive application due to their low water solubility and difficulties of their removal from soil. To increase the water solubility of the cyclodextrins and their complexes, Petrovi et al. investigated many substituted derivatives prepared by the reaction of β-cyclodextrin with different reagents. Methylation of β-CD increases its' aqueous solubility, extends the hydrophobic cavity, and provides a greater surface area for the complexion of pesticides like dimethoate. The water-soluble methylated-β-cyclodextrin may be potentially useful for improving the application of hydrophobic organic pesticides and enhancing their removal from the environment and remediation of contaminated soil and groundwater [34].
Furthermore, the authors omitted to clarify the mechanism of inclusion with β-CD despite the variety of characterization used approaches gain depth in sight into the inclusion of complex formation between dimethoate and β-CD. Until now, no theoretical investigation has used DFT calculations to evaluate the interactions between dimethoate and β-cyclodextrin and their nature.
Thus, this study is provided in gas and aqueous phases to describe the nature of non-covalent inter-molecular interactions, namely, the H-bond ones between DMT and β-CD, during the inclusion of complex formation.
To identify the nature of the interaction between the host and the guest, we applied natural bond orbitals (NBO), energy decomposition analysis (EDA), quantum theory of atoms in molecules (QTAIM), non-covalent interactions-reduced density gradient (NCI-RDG), and independent gradient model (IGM).

Computational methods
DFT calculations were carried out with the Gaussian 09 software package [35], and the creation of molecular graphs was carried out by the GaussView 6 programs [36].
The structure of dimethoate was obtained for the first part of our investigation from the PubChem compound database [37]. Additionally, ChemOffice 3D Ultra was used to build the structure of the β-CD (version 10, Cambridge software) [38]. After that, each configuration was fully optimized using the MOPAC 2016 package's semi-empirical approach PM7 [39].
The initial inclusion structures of DMT and β-CD were constructed by the HyperChem version 7.5 molecular modeling program [40]. The guest dimethoate was moved to enter and cross the host cavity of β-CD along the Z-axis from − 6 to 6 with one pass in each of the two possible directions, A and B (see Fig. 2). The PM7 approach was used to optimize the created structures at each step without any restrictions, and the most stable complexes were identified. The created structures at each step were optimized at the PM7 method without any symmetrical restrictions deducing the most stable complexes.
The solvent effect on conformational equilibrium was introduced in the single-point DFT calculation of water (ɛ = 78.5) using the conductor polarizable continuum model (CPCM) [50]. Global reactivity indices were determined using the HOMO-LUMO transitions reactivity indices obtained through HOMO-LUMO transitions.
According to the minimum energy structure, the complexation energy upon complexes between dimethoate and beta-cyclodextrin is given by Eq. (1) [51].
where, E complex , E freeDMT , and E freeβ-CD correspond to the complex's heat formation (HF) energies, the free guest DMT and the free host β-CD, respectively.
The intensity of the energy variation would be a sign of the force driving toward complexes.
The interaction and deformation energy (DEF) for each host and guest component during complex formation is the difference between the optimized component compared to its energy in the complex (Eq. (2)) [52]. where E (component)spopt is the single-point energy of the component using its geometry in the optimized complex and E (component)opt is the energy of the optimized geometry of the component. The difference between the energy of the fully optimized component and its energy in the complex was used to determine the deformation energy for each component, host, and guest during the complex's formation [53].
QTAIM, NCI-RDG, and IGM analysis were explored using the Multiwfn program [58] and visualized by the VMD program [59].

Energetic and structural analysis of complexes
The PM7 calculation was established for models A and B to control the complexation energy during the inclusion process of DMT in the β-CD cavity from − 7 to + 7 to identify global minimum energy. Figure 2 shows the variation of the complexation energy (ΔE) during the inclusion process of both models, as A and B are a function of the Z distance. These results indicate that the complexation energy values are negative, implying that the DMT@ β-CD-produced complexes are energetically favorable [59,60]. For both A and B orientations, the global minimums of the most stable structures are situated at Z = − 2 Å. From Fig. 3, it can be seen that the complexation energy values obtained are − 42.86 kcal/mol and − 44 kcal/mol for orientations A and B, respectively.
PM7 calculations evaluating the complexation energies show that orientation B gives rise to the conformer, which is more stable by 1.14 kcal/mol than that corresponding to orientation A.
The calculated complexation and interaction energies in gas and aqueous phases with the functional CAM-B3LYP for orientations A and B are provided in Table 1.
However, orientation B provides the most significant values of E encapsulation. Aqueous phase interaction  energies lend support to gas phase interaction energies. Compared to orientation A, orientation B has stronger negative energy. Additionally, we can see that the interaction energy difference in the gas and aqueous phases, respectively, is 1.48 kcal/mol and 8.708 kcal/mol, supporting the stability of orientation B.
The results also revealed that the ΔE interaction energy of the β-CD molecule with functional CAM-B3LYP/6-31G (d, p) is higher than that of the DMT molecule in the two orientations A and B. E interaction is a crucial measure in assessing the stability of inclusion complexes.
In addition, ΔE interaction is an important parameter measuring the stability of inclusion complexes; the results reveal that the deformation energy of the β-CD molecule with functional CAM-B3LYP/6-31G (d, p) is higher than that of the DMT molecule for both orientations. This interaction energy demonstrates that the β-CD structure's flexibility is crucial for increasing intermolecular interaction and the entire system's stability after complexation. Figure 4 shows the favorable structures of the lowest energy conformers obtained in the water and gas phase using the CAM-B3LYP functional. It is clear from the structures that the DMT molecule is completely encapsulated in the β-CD cavity when the intermolecular hydrogen bonds (HBs) in both structures control the stabilization between the two orientations, A and B, in different phases.

TD-DFT analysis
The UV-vis spectra are obtained using the DT-DFT method with CAM-B3LYP functional and 6-31G (d,p) basis set in water to identify variation electron transitions [61]. The absorption wavelength, oscillator strength (f), minor and major orbital contributions, and their predicted energies (E) have been reproduced and shown in Fig. 5 and Fig. 6. The number of transitions for calculation for absorption spectra was 50 states.
The HOMO and LUMO orbitals of the two orientations are entirely localized on the DMT. These results revealed that encapsulation does not change the guest molecule's charge distribution. Analysis of the simulated absorption spectra of the inclusion complexes revealed that for orientation A, the maximum of the absorption λ abs is 172.15 nm, which  corresponds to the (HOMO − 1) → (LUMO + 1), which was ascribed to LE(DMT) with a contribution of 45.96% and oscillator strength f = 0.075. Moreover, for orientation B, the maximum simulated peak was discovered theoretically at 180.67 nm due to the contributions of 30.7% (HOMO − 9) → (LUMO), which were signed to CT (β-CD → DMT). Both orientations' UV-vis absorbance after complexation with β-CD differs from free DMT. It has been proposed that these variations in the absorption bands are generated by conformational changes in the structure of DMT that occur when the inclusion complex with β-CD is formed.
Both orientations' UV-vis absorbance changes after complexation with β-CD. Variations in absorption bands may be due to conformational changes in DMT during complex inclusion formation with β-CD.

Molecular reactivity analyses
The HOMO and LUMO are crucial chemical indicators largely used in this investigation. The maximum electronic charge (ΔN), electronic potential (μ), hardness (η), and global electrophilicity index (ω) as global indices of reactivity are calculated using the following equations [62][63][64]: where ΔN(MAX)host or guest = (− µ/η) host or guest The computed HOMO and LUMO energies and reactivity parameters of inclusion complexes in the gas phase and water are calculated at CAM-B3LYP /6-G (d, p) are gathered in Table 3.
From Table 3, we noticed that the HOMO-LUMO energy gap value in the water phase is 8.73 eV for orientation A and 9.32 eV for orientation B. The value of ( E HOMO − E LUMO ) gap for orientation B is higher than that of orientation A, which indicates that orientation B gives rise to a more stable complex than that of orientation A. This result and the calculated binding energy are in good agreement.
The chemical potential of both orientations is negative, leading to a spontaneous inclusion process is spontaneous. μ free gues t ˃ μ free host ; this indicates that the direction of the charge transfer associated with the creation of the inclusion complex is from DMT to β-CD. We observed that the most significant chemical hardness (η) value for model A is 4.62 eV in water, comparable to that of orientation B (4.49 eV), indicating that the charge transfer in orientation B is significant. Compared to orientation A, orientation B exhibits a higher global electrophilicity index (ECT) value, indicating that charge transfer happens from the host to the guest. The result shows that orientation B corresponds to the most electrophilic complex.
For each pair i (donor) and j (acceptor), the stabilization energy E (2) accompanying delocalization i → j is well defined in the literature [71,72], which is given by Eq. (9).
qi is the donor orbital occupancy, εi and εj are diagonal elements, and F(i,j) is the off-diagonal NBO Fock matrix elements.
The stabilization energy E (2) and bond length related to the most considerable interactions for complexes B identified by the CAM-B3LYP/6-31G (d, p) method for B orientation in both gas and water phase are summarized in Table 4.
According to the results grouped in Table 4, two different classes of categories exist weak hydrogen formed between (LP) and (BD*) and van der Waals interaction created between the (BD) and (BD*).
The higher stabilization second-order perturbation energy E (2) is 1.28 kcal/mol, associated with high donor-acceptor interaction corresponding to an hydrogen bond length of 1.83 Å.
Accordingly, NBO calculations highlight the hydrogen bond's contribution to sustaining the host-guest reaction and maintaining stability.

Quantum theory of atoms in molecules (QTAIM)
The QTAIM analysis plays a crucial role in identifying intra-and intermolecular interactions. Quantum mechanical parameters such as electron density at the critical bond points (BBCPs) are used to determine the nature of host-guest interactions and classify bonding interactions [73].
The main topological parameters to define the properties of critical bond point BCPs are the total electron density ρ(r) and its Laplacian ∇ 2 ρ(r) [74][75][76].
Thus, Table 5 collects the QTAIM characterizing parameters of the (3, − 1) critical points of the DMT@β-CD complex. The QTAIM molecular graphs representing orientations B are illustrated in Fig. 7.
The results gathered in Table 5 of the QTAIM calculation of model B using the CAM-B3LYP/6-31G (d, p) functional in gas and water shows an interaction between dimethoate and β-CD through an H-bond.
From CAM-B3LYP/6-31G(d, p) results, ρ(r) values vary from 0.001 to 0.03 a.u and 0.001 to 0.029 a.u, respectively, for gas and water, while Laplacian ∇ 2 ρ(r) values are in the range from 0.004 to 0.094 in gas and 0.0059 to 0.096 a.u in water.
The results of ρ(r) in the gas phase and water showed values in the range of 0.001 to 0.03 a.u and of 0.001 to 0.029 a.u, respectively, with the corresponding Laplacian ∇ 2 ρ(r) varying between 0.004 and 0.094 a.u in the gas phase and 0.0059 and 0.096 a.u in water.
However, stronger hydrogen bonding O153--H144 is observed with the lowest intermolecular distance of 1.83 Ǻ in the gas and aqueous phase and the maximum electron density ρ(r) and Laplacian ∇ 2 ρ(r). The ellipticity values for the intermolecular bonding of the DMT@β-CD complex range from 0.002 to 1.23 a.u in gas and from 0.02 to 0.23 a.u in water, indicating stable contact between the host and guest [81].
All calculated ∇ 2 ρ(r) and H(r) values are positive, indicating the presence of weak electrostatic interactions, and the calculated ratio of − G(r)/V(r) is > 1 relative to significant interactions of the non-covalent character.
The topological parameters Table 5 in the gas phase and water display all H(r) values are positive, and ∇ 2 ρ(r) are all small positive values implying weak interaction mainly of electrostatic. Besides, the ratio of − G(r)/V(r) is > 1 for the complex, supporting the existence of weak intermolecular bonding. The Table 5 bond energy (E) values show that the principal molecular interaction in the gas phase or water is detected for O153--H144 with − 0.011 kcal/mol.
The result of ʎ1, ʎ2, and ʎ3, corresponding to the Hessian eigenvalues of the electron density at BCP, indicates that ʎ1 ˂ ʎ2 ˂ ʎ3, the sum of negative curvatures (ʎ1 + ʎ2) as well as the positive. At the same time, (ʎ3) decreases with H⋯O distances.
This result shows that the electron density increase in the plane perpendicular to the bond path occurs concurrently with electron density depletion along the bond path. Lower ellipticity index values demonstrate that electrons are delocalized through the associated atoms.
The QTAIM results show that van der Waals interactions and weak hydrogen bonds are the chief factors influencing the complex's stability.

NCI-RDG analysis
To identify hydrogen bonds, van der Waals contacts, and repulsive steric interactions between host and guest in the formed complex, the non-covalent interaction (NCI) via a reduced density gradient (RDG) was employed [82,83].
The equation that describes the RDG approach is as follows [82]: Considering Fig. 8a, which plots RDG against sign (λ 2 ) ρ; the sign (λ 2 ) value can be used to indicate the type of interaction; for example, sign (λ 2 ) > 0 indicates a repulsive interaction, while sign (λ 2 ) < 0 indicates an attractive interaction, such as hydrogen bonds.
It is observed that van der Waals interactions are ranged from − 0.018 to 0.005 a.u and are shown with a green spot; the hydrogen bonding interactions are illustrated with a blue spot and located between − 0.05 and − 0.02 a.u. The red spot indicates the repulsive steric forces.
From the 3D spatial NCI isosurface diagram (Fig. 8b), we can see critical green patches in the region between DMT and β-CD related to van der Waals interactions, indicating that the guest forms a stable inclusion complex with the host. Besides, a blue patch represents the strong hydrogen bonding interactions and red spots represent the repulsive steric forces.
The resulting 2D scatter plot shows that red points correspond to δ ginter , while black points represent δ gintra . Figure 9, with the most intense black peak appearing on the negative side at sign (λ 2 )ρ = − 0.28 with δ ginter of 0.392 au. van der Waals interactions can be seen in the region where sign (λ 2 )ρ = − 0.04, with δ ginter of approximately 0.056 a.u. corresponding at the second less intense peak. In contrast, the positive side of the sign (λ 2 )ρ peak lies in the range of 0.04 to 2.00 with δ ginter = 0.168 a.u., indicating a repulsive interaction.
The 3D IGM isosurface map for encapsulated complexes is depicted in Fig. 9. The green-colored regions represent weak van der Waals interactions, whereas blue regions denote stronger electrostatic attraction. The results show that the DMT@β-CD complex is stabilized by hydrogen bonds and van der Waals interactions. From these results, it can be concluded that there are intermolecular hydrogen bonds, and the IGM analysis agrees with the QTAIM results.

Energy decomposition analysis (EDA)
To highlight and evaluate the hydrogen bonding between DMT and β-CD in gas and aqueous phases, the Morokuma and Ziegler-Rauk energy decomposition analysis (EDA) [86][87][88] was applied, which was widely used previously [89][90][91][92][93][94][95][96]. Thus, the EDA using the hybrid B3LYP-D3 functional gives rise to the resulting interaction energy ΔE int which is decomposed into four terms of energy as given below: the electrostatic interaction ΔE elstat which is an attractive interaction, the Pauli interaction ΔE Pauli that exhibits a repulsive interaction, the orbital interaction ΔE orb as an attractive term refers to the charge transfer between the occupied orbitals, and the unoccupied orbitals of the two fragments and finally ΔE disp corresponding to the Grimme dispersion correction term.
(10) ΔE int = ΔE elstat + ΔE Pauli + ΔE orb+ ΔE Disp  The ∆E int energies gathered in Table 6 and Fig. 10 are negative, indicating the stabilization effects of these interactions. The analysis of the results shows that the interaction energy of − 43.88 kcal/mol in the gas phase is higher than that obtained in the water phase of − 70.2 kcal/mol. The divergence between the two values comes essentially from the ΔE Pauli , ΔE elstat , and ΔE orb contributions to the ∆E int total interaction energy, where the former is considerably weakened. In contrast, the two later terms are substantially strengthened (increased absolute values), as summarized in Table 6. Indeed, the ΔE orb goes up from 20 to 35% and increases from 4 to 44%. Besides, the changing middle − 34.73 vs. − 33.88 kcal/ mol does not influence the contribution. The ΔE disp relative to the hydrogen interactions contributes 39.12% in gas and 31.77% in water into the total interaction energy. Thus, emphasize the strong hydrogen interaction type between DMT and β-CD.

Conclusion
This theoretical reports the interactions between DMT and β-CD in gas and aqueous phases. The results indicated that orientation B is more privileged than A in both cases, where the optimized structures of the DMT@β-CD complex confirm the total inclusion of dimethoate into the β-CD cavity. NBO results revealed that the guest molecule interacts with the host and modifies its charge distribution to form a stable inclusion complex. QTAIM, NBO, and EDA analyses confirm the existence of hydrogen bonds between β-cyclodextrin and dimethoate. RDG and IGM calculations showed that DMT@β-CD forms a combination of weak hydrogen bonds and van der Waals interactions stabilizing the complex. NCI isosurface confirms the establishment of hydrogen bonds and, vdW, and steric and steric repulsion during complex inclusion formation. Considering the inclusion of dimethoate in β-cyclodextrin, this study could serve as a starting point for experiments on pesticide environmental problems. The UV-vis electronic spectra shows significant electronic transitions giving rise to inter-ligands charge transfers (ILCT).