Integrating molecular modeling methods to study the interaction between Azinphos-methyl and gold nanomaterials for environmental applications

1.   Introduction

Azinphos-methyl (AzM) is an organophosphorus compound with the molecular formula (C₁₀H₁₂N₃O₃PS₂). It has a density 1.44 g/cm³ and a molar mass 317.32 g/mol ^[1]. AzM is a pesticide that is commonly used in agriculture to control a variety of pests, including fruit flies, moths, and caterpillars. It belongs to a class of pesticides known as organophosphates, which operate by inhibiting the activity of an enzyme called acetylcholinesterase, leading to an accumulation of acetylcholine in the nervous system, causing nerve damage and ultimately resulting in the death of the pest ^[2]. Investigating the interaction among metal nanoparticles and molecules is an important area of research in materials science and nanotechnology ^[3]. The motivation for studying the interaction between Azinphos-methyl (AzM) and metal nanoparticles, such as gold nanoclusters (AuNCs) and nanoparticles (AuNPs), arises from the need to enhance detection and analysis of pesticides using advanced nanomaterials. Metal nanoparticles exhibit unique optical properties, like localized surface plasmon resonance, which significantly enhance Raman signals in surface-enhanced Raman spectroscopy (SERS) that make them attractive for a wide range of applications, including sensing, catalysis, and biomedicine ^[4]. This enhancement allows for sensitive detection and detailed characterization of AzM, improving our understanding of its behavior and interactions. However, to fully realize the potential of metal nanoparticles, it is important to understand how they interact with molecules that can affect their behavior and properties ^[5]. The interaction between metal nanoparticles and molecules can be studied using various techniques, including spectroscopy, microscopy, and electrochemistry. For example, surface-enhanced Raman spectroscopy (SERS) ^[6] is a powerful technique that can be used to study the interaction between molecules and metal nanoparticles ^[7]. SERS involves adsorbing molecules onto the surface of metal nanoparticles and measuring the enhancement of Raman scattering due to the interaction between the molecule and the nanoparticle's surface ^[8]. Other techniques, such as transmission electron microscopy (TEM) and atomic force microscopy (AFM), can be used to directly visualize the interaction between metal nanoparticles and molecules ^[9]. The structural and electronic properties of AzM containing gold clusters were investigated using density functional theory (DFT)/B3LYP with 6-311++G/LANL2DZ mixed basis sets ^[10]. Both atoms-in molecules theory (AIM) and Natural bond orbital (NBO) analysis have been applied to a wide range of chemical systems, including small molecules, proteins, and materials. These methods are useful tools for predicting the properties and reactivity of molecules, as well as for understanding the fundamental principles that govern chemical bonding and reactivity ^[11]. The interactions of AzM with AuNPs of different sizes and shapes were investigated using Monte Carlo (MC) simulation ^[12]. Numerous theoretical studies have explored the interactions between molecules and gold nanoclusters (AuNCs) or nanoparticles (AuNPs). For instance, Mageed and Taha ^[13] investigated how serine and threonine dipeptides interact with gold nanoclusters and nanoparticles of varying shapes and sizes. Their research used quantum mechanics and molecular simulations to quantitatively assess adsorption energies, revealing significant insights into peptide adsorption preferences and interactions in different phases (gas and liquid). These findings are crucial for understanding peptide and protein behavior in the presence of gold nanomaterials, with implications for nanotechnology applications. In another study, Shweta Bhardwaj and colleagues ^[14] analyzed interactions between molecules and small gold clusters (Au_n, n = 6–9) using density functional theory (DFT). They calculated geometric parameters, interaction energies, and performed natural bond orbital (NBO) analysis to examine charge transfer and molecular orientation effects on Raman enhancement. Their results quantitatively demonstrated how variations in gold cluster size impact Raman activity and stability. Hariharan and Vadlamudi conducted three separate studies ^[15,16,17] on the interaction between molecules and small gold clusters (Au_n, n = 1–4). Their work quantitatively showed that the number of atoms in the gold cluster correlates with enhanced Raman activity and stability, and they evaluted interaction energies, reactivity parameters, electrostatic properties, and used NBO analysis. Silver clusters were utilized ^[18] for the detection of the pterin. Our aim is to identify the interactions between Ptr and silver NCs. In particular, the binding energies between Ptr and Agnq (n = 1–6; q = 0, +1, +2) NCs, established the complexes that are the most likely to be formed in aqueous solutions, analyzed the nature of the bonds between silver and Ptr, calculated absorption and Raman spectra of silver-Ptr complexes. This work ^[19] surface enhanced Raman spectroscopy (SERS) experiments, and density functional theory (DFT) calculations were carried out to investigate adsorption of cysteine on gold nanoparticle surface. Therefore, in this work, the conformation and adsorption sites of the flexible typical amino acid cysteine on gold nanoparticle surface were determined by combining both DFT and SERS methods.

In recent decades, SERS spectroscopy has proven to be fundamentally important in various research fields, including chemical sensors, materials science, and biomedicine. For these applications, detailed information on molecule-metal interactions is essential ^[20]. The combined SERS-DFT-MC approach used in this study provides a comprehensive view of molecular interactions with metal surfaces.

2.   Computational details

2.1.   DFT calculations

DFT calculations for geometry optimization with two alternative basis sets were performed using the Gaussian 09 program ^[21]. The 6-311++G (d, p) basis set was employed for C, H, N, P, S, and O atoms, while the double-z basis set of Hay and Wadt (LANL2DZ) was used for Au atoms. The Becke's three-parameter Lee-Yang-Parr exchange-correlation method was utilized for the optimization process ^[22]. The functional used for investigating the electronic properties and interactions between Azinphos-methyl and gold clusters was the B3LYP functional. This functional was suitable for conducting vibrational spectroscopic analysis and generated theoretical results that align with experimental observations ^[23]. Optimized geometries were identified as minima on the potential energy surface, and all vibrational frequencies exhibited positive values. Individual optimizations were first conducted for both the molecules and Au clusters. Subsequently, optimizations were performed with the molecules present in the presence of the gold clusters ^[24]. The Raman activity was computed based on the optimized geometries. The Raman intensities, which were initially calculated using Gaussian 09, were later converted into relative Raman activities. Finally, the spectral data were visualized using the Origin software ^[25].

By utilizing the HOMO (highest occupied molecular orbital) and LUMO (lowest unoccupied molecular orbital) energy values of a molecule, along with parameters such as ionization potential ${I}_{P}$ and electron affinity ${E}_{A}$, electronegativity (χ), softness (σ), hardness (η), and global electrophilicity (ω), the calculations can be performed as outlined below ^{[21,22,23,24]}:

Ionization potential (${I}_{p}$) and electron affinity $\left({E}_{A}\right)$ are two related reactivity parameters that describe the tendency of an atom or molecule to lose or gain electrons, respectively.

Electronegativity is a measure of the tendency of an atom to attract electrons towards itself when bonded to another atom. It is a fundamental property of atoms that influences their chemical behavior and reactivity.

Hardness typically refers to the resistance of a material to deformation or indentation. It is expressed as:

Inverse of hardness is the global softness, ${\rm{ \mathsf{ σ} }}$, defined as:

Hardness and softness are related to the ability of a molecule to undergo chemical reactions. Hard molecules have high ionization energies and low electron affinities, making them less reactive and more stable, and have a large energy gap. Soft molecules, on the other hand, have low ionization energies and high electron affinities, making them more reactive and less stable, and have a small energy gap ^[26].

The electrophilicity index is given by:

The electrophilicity index is a measure of the electrophilicity of a molecule or a part of a molecule, which is related to its ability to accept electrons in a chemical reaction. It is a quantitative measure of the tendency of a molecule to act as an electrophile and is calculated based on the chemical potential and hardness of the molecule ^[27].

The dipole moment µ is a measure of the separation of positive and negative charges within a molecule. It quantifies the overall polarity of a molecule and is represented by a vector pointing from the negative to the positive charge ^[28].

The static polarizability (α₀) is a measure of how easily a molecule can be polarized by an external electric field. It describes the ability of a molecule to undergo a distortion in its electronic charge distribution in response to an applied electric field ^[29].

The total polarizability (Δα) for each complex was calculated using the [14, Equation 8]:

The concept of natural charge transfer (Δq) ^[30] is important in understanding the electronic properties of molecules and their behavior in chemical reactions and interactions. It helps to describe the flow of electrons within a molecule and can provide insights into the reactivity, stability, and properties of molecules $\varDelta q$, defined as.

The molecular electrostatic potential (MEP) is a representation of the electrostatic potential energy experienced by a positive test charge placed in the vicinity of a molecule ^[31,32]. It provides information about the distribution of positive and negative electrostatic potentials around the molecule and can help visualize the regions of high and low electron density ^[15].

The binding energy for each complex was calculated using the Eq 10:

where ${E}_{({Au}_{n}-AzM)}$, ${E}_{{Au}_{n}}$, and ${E}_{AzM}$ are the total energy of the complex, AuNCs, and AzM, respectively, the binding energy represents the strength of the interaction between the molecule and the gold cluster. A negative binding energy indicates that the interaction is favorable, meaning that the molecule and the gold cluster are more stable when they are bound together than when they are isolated. A positive binding energy indicates that the interaction is unfavorable ^[33].

NBO analysis of the charge population was conducted for the optimized structures to investigate the charge delocalization between the AuNCs and AzM by utilizing second-order perturbation theory. The stabilization energy E² related to electron delocalization between an acceptor orbital and a donor orbital is an important concept in NBO analysis. It is determined as

where ${q}_{i}$ is the occupancy of the ith donor orbital. The term ${F(i.j)}^{2}$ refers to the off-diagonal elements of the kohn-Sham matrix in the NBO analysis. The ${\varepsilon }_{i}$ and ${\varepsilon }_{j}$ are the energies of the donor and acceptor orbitals ^[34], respectively.

The AIM theory ^[11] was created to improve our understanding of chemical bonding by examining the electron distribution within AzM-Au_n complexes. This was done through a detailed topological examination at the bond critical point (BCP) using the Multiwfn software ^[35], based on the wave function obtained from density functional theory (DFT) calculations. The criteria for evaluating pre-existing associations included examining ρ(r) and ∇²ρ(r) at the BCP. The (H(r)) electronic energy density at the BCP is determined by combining the (G(r)) kinetic energy and (V(r)) potential energy ^[13].

2.2.   Monte Carlo simulation

MC simulation is a computational technique that can be utilized to simulate the adsorption of a molecule onto a surface and calculate the adsorption energy ${(\varDelta E}_{ads})$ ^[36].

The adsorption energies of AzM on a simple box, right circular cone, cylinder, and sphere gold nanoparticles at different sizes (3, 5, and 7 nm) were investigated by the Materials Studio (BIOVIA, 2020) package ^[37]. First, the materials studio database was used to extract a cubic gold crystal. The gold unit cell was a face-centered cubic (FCC) structure with a length of 4,078 Å. The FCC structure's space group symmetry for cubic gold crystals was identified as Fm-3m (225), belonging to the m-3m crystal class. The structure underwent optimization using the COMPASS Ⅱ force field ^[38], a method implemented through the Forcite Module ^[39]. The convergence tolerance for energy, stress, displacement, and force were 2.e^-5 kcal mol^-1, 0.001 GPa, 1.0e^-5 Å, and 0.001 kcal mol^-1 Å^-1, respectively. The AuNPs shapes were made using build = > build Nanostructure = > Nanocluster tool, as follow: Cubic box, 3nm (X, Y, Z directions = 30 Å), 5 nm (X, Y, Z directions = 50 Å), 7 nm (X, Y, Z directions = 70 Å), Right circular cone, 3 nm (Radius = 15 Å, Height = 30, Base plane (1 1 1)), 5 nm (Radius = 25 Å, Height = 50, Base plane (1 1 1)), 7 nm (Radius = 35 Å, Height = 70, Base plane (1 1 1)), Cylinder, 3 nm (Radius = 10 Å, Height = 30, Base plane (1 1 1)), 5 nm (Radius = 10 Å, Height = 50, Base plane (1 1 1)), 7 nm (Radius = 10 Å, Height = 70, Base plane (1 1 1)), Sphere, 3 nm (Radius = 15 Å), 5 nm (Radius = 15 Å), and 7 nm (Radius = 15 Å). The AzM was also optimized by the same Module that was used in the gold structure ^[13]. The forcefield COMPASS Ⅱ (Condensed-phase optimized molecular potentials for atomistic simulation studies) was assigned for the MC simulation, which was conducted by the Adsorption Locator Module in the BIOVIA Materials Studio 2020 package ^[40]. This force field is well-suited for studying adsorption processes, which are crucial in various fields like catalysis, materials science, and environmental science. It incorporates parameters for non-bonded interactions within the molecule and between the molecule and the surface, enabling the calculation of potential energy during adsorption. The disparity in potential energy between the adsorbed and unadsorbed states is then used to determine the adsorption energy ^[37]. The Monte Carlo simulation uses simulated annealing with 10 cycles, each with 100,000 steps. The Metropolis Monte Carlo technique ^[41] is utilized by the Adsorption Locator Module to determine the adsorption configuration with the lowest energy.

3.   Results and discussion

3.1.   DFT calculations

The properties of a molecule are intrinsically linked to its optimized structure, which represents the configuration with the minimum energy. Identifying this relaxed structure is crucial as it determines the molecule's stability and reactivity. Figure 1 displays the optimized structure of Azinphos-methyl (AzM), obtained through geometry optimization using the B3LYP method with a 6-311++G (d, p) basis set. The B3LYP method, a widely used density functional theory approach, provides accurate predictions of molecular structures by minimizing the total energy of the system. The 6-311++G (d, p) basis set, with its diffuse and polarization functions, ensures a precise description of electron distribution, a crucial for evaluating molecular interactions and properties ^[22].

Figure 1.  DFT-B3LYP/6-311++G optimized structure of AzM.

DownLoad: Full-Size Img PowerPoint

3.1.1.   Comparison of experimentally and computationally SERS spectra of AzM

In this section, we present the Raman spectrum calculated using density functional theory (DFT) for AzM. Figure 2 shows a notable resemblance between the calculated spectrum and the experimental data reported by Liu et al. ^[42]. We have also provided a comparison of the Raman peaks derived from DFT calculations with those observed experimentally in Table 1. Moreover, the table shows the vibrational assignments for these peaks, which were calculated by direct observation using GaussView6 software ^[43].

Figure 2.  DFT calculated the Raman spectrum of AzM.

DownLoad: Full-Size Img PowerPoint

Table 1.  The comparison between the experimental Raman peaks of AzM and the DFT calculated Raman peaks of AzM, and vibrational mode Assignment for AzM.

No	Raman peaks		Assignment of Raman peaks
	Computationally	Experimental [42]
1	593	587	δ(C=O) deformation vibration
2	683	674	υ(P=S) stretching
3	714	703	Benzene ring breathing
4	791,868	775,897	1,2,3-triazine ring breathing
5	1057	1026	Asymmetric P-O-C deformation vibration
6	1127	1221	γ(C-H) in P-O-CH3
7	1239	1258	υ(C-N) in S-CH2-N
8	1288	1283	υ(C-N) in O=C-N
9	1309	1332	1,2,3-triazine ring breathing
10	1365	1399	γ(C-H) in S-CH2-N
11	1453	1450	υ(N=N) stretching
12	1491	1495	1,2,3-triazine ring breathing

---

 | Show Table

DownLoad: CSV

According to the research, there is a significant similarity between the peaks obtained through density functional theory (DFT) and those seen in the Surface-Enhanced Raman Scattering (SERS) spectrum. However, there are slight differences or inconsistencies between the DFT-calculated peaks and those obtained from the SERS spectrum ^[42]. These differences can be attributed to an acceptable margin of error, which is illustrated in Table 1.

3.1.2.   Optimization of Au_n (n=2, 4, and 6) clusters using DFT calculations

To optimize Au_n clusters (n = 2, 4, and 6), we used density functional theory (DFT) calculations. The aim is to find the most energetically stable arrangements of gold clusters that contain 2, 4, and 6 atoms, using the B3LYP/LANL2DZ level of theory. The optimization process involved adjusting the positions of the atoms within the cluster to determine the minimum energy configuration. This adjustment was done by solving the Schrödinger equation using DFT, which provides a theoretical prediction of the electronic and geometric structure of the cluster. Once the optimization process is complete, we can analyze the optimized geometries to investigate the electronic and geometric properties of the different gold clusters, as well as their potential interactions with other molecules or surfaces. The optimized geometries of the different gold clusters are presented in Figure 3.

Figure 3.  Optimized Geometries of Au_n (n = 2, 4, and 6).

DownLoad: Full-Size Img PowerPoint

3.1.3.   Geometric analysis of AzM near Au_n (n=2, 4, and 6) clusters using DFT: Optimization and Raman activity

To optimize AzM-Au_n (n = 2, 4, and 6), we used density functional theory (DFT) calculations using the three-parameter hybrid-functional of Becker (B3LYP), with a 6-311++G (d, p) basis set was used for C, H, N, P, S and O atoms and LANL2DZ for the Au atoms (Figure 4).

Figure 4.  The optimized geometries of AzM-Au₂, AzM-Au₄, and AzM-Au_6, respectively.

DownLoad: Full-Size Img PowerPoint

In Figure 5, a comparison is made between the calculated Raman spectra of the AzM molecule and AzM-Au_n (n = 2, 4, and 6) complexes. The spectra exhibit five distinct peaks, demonstrating the anticipated pattern of increased Raman activity as the number of Au atoms in the cluster rises. These peaks, characterized by average wavenumbers around ~594, ~683, ~867, ~1239, and ~1453 cm⁻¹ are detailed in Table 2, along with their assigned vibrational modes. The enhancement factors, defined as the ratio of the Raman intensity of AzM-Au_n to that of AzM, follow the same trend as the Raman activities, indicating a correlation between the number of gold atoms and the increased Raman signal. Specifically, the vibrational band at 594 cm⁻¹ corresponds to the δ(C = O) deformation vibration. Moreover, the bands at 683 cm⁻¹ and 867 cm⁻¹ are linked to the υ(P = S) stretching and 1,2,3-triazine ring breathing, respectively. Additionally, the bands at 1239 cm⁻¹ and 1453 cm⁻¹ correspond to the υ(C-N) in S-CH₂-N and υ(N = N) stretching, respectively. AzM-Au₆, however, does not exhibit the anticipated rise in Raman activity, which will be further discussed below.

Figure 5.  The calculated spectra of (a) AzM, (b) AzM-Au₂, (c) AzM-Au₄, and (d) AzM-Au₆.

DownLoad: Full-Size Img PowerPoint

Table 2.  The enhancement factors (EF) of Raman activity (RA) for different peaks of the AzM-Au_n complexes, along with their corresponding vibrational mode assignments.

Peak (cm-1)	AzM RA	AzM-Au2		AzM-Au4		AzM-Au6		Assignments
		RA	EF	RA	EF	RA	EF
594	11.59	50.50	4.35	222.96	19.23	133.89	11.55	δ(C=O) deformation vibration
683	48.07	174.99	3.64	278.83	5.80	257.07	5.34	υ(P=S) stretching
867	20.24	50.75	2.50	93.59	4.62	81.52	4.02	1,2,3- triazine ring breathing
1239	96.58	120.92	1.25	170.28	1.76	144.46	1.49	υ(C-N) in S-CH2-N
1453	59.73	83.17	1.39	114.52	1.91	107.27	1.28	υ(N=N) stretching

---

 | Show Table

DownLoad: CSV

3.1.4.   The electronic properties

The electronic properties of AzM and AzM-Au_n (n = 2, 4, and 6) were calculated by the DFT-B3LYP with 6-311++G (d, p) and LANL2DZ basis set using the Gaussian 09. The analysis encompassed a range of electronic characteristics, including E_HOMO, E_LUMO, ΔE_gap, I_P, E_A, χ, $\eta$, $\mathcal{w}$, $\sigma$, ΔE_bind.

When scientists want to understand how a molecule will act in a chemical reaction, they look at two special parts of the molecule called the E_HOMO and E_LUMO. HOMO (highest occupied molecular orbital) refers to the electron-donating (nucleophilic) nature of a molecule. The higher the energy of the HOMO, the easier it is for the molecule to release electrons. LUMO (lowest unoccupied molecular orbital) refers to the electro-acceptor (electrophilic) character of a molecule. The lower the energy of the LUMO, the easier it is for the molecule to accept electrons. Figure 6 illustrates the electronic structures of the HOMO and LUMO for pure AzM and AzM-Au_n (n = 2, 4, and 6).

Figure 6.  Visualization of the HOMO and LUMO levels of (a, b) AzM, (c, d) AzM-Au₂, (e, f) AzM-Au₄ and (g, h) AzM-Au₆.

DownLoad: Full-Size Img PowerPoint

The energy gap (E_g) between the HOMO and LUMO energy levels in a molecule is another parameter that can be used to predict the degree of stability of a chemical entity. High values of the energy gap relatively imply high stability of the molecule in a given environment ^[22], and lower values indicate that the molecule is highly active in a chemical reaction. The most stable chemical compound is the one with the greatest energy difference between HOMO and LUMO. By looking at these parts of the molecule, scientists can determine how it will react with other substances ^[44].

The results in Table 3 indicate that the gap energy of the AzM-Au₄ complex is the lowest compared to the other complexes, namely AzM-Au₂ and AzM-Au₆, and that it will therefore be the least stable. Second, these results indicate that the AzM-Au₂ and AzM-Au₆ complexes have similar stabilities, requiring further approaches and calculations to determine which of these complexes is more stable. In addition, all complexes show a lower energy gap than the AzM molecule, an overall property in which the new molecules are more flexible and possess novel electronic properties.

Table 3.  The E_HOMO and E_LUMO and E_g for the AzM and AzM-Au_n (n = 2, 4, and 6).

	Electronic states		Eg (eV)
	EHOMO (eV)	ELUMO (eV)
AzM	–6.77	–2.60	4.16
AzM-Au2	–5.70	–2.93	2.77
AzM-Au4	–5.41	–3.62	1.79
AzM-Au6	–5.87	–2.97	2.89

---

 | Show Table

DownLoad: CSV

The data presented in Table 4 were obtained using Koopman's theorem, which calculates a molecule's ionization potential and electron affinity using its HOMO and LUMO energy levels. The ionization potential measures the energy required to remove an electron from a molecule, while the electron affinity measures the energy released when adding an electron to a neutral molecule, resulting in the formation of a negative ion. Table 5 provides various reactivity parameters, such as Hardness, Softness, and Electrophilicity. The global electrophilicity index ($\mathcal{w}$) values show an increasing trend as the number of atoms in the gold cluster increases. For instance, in AzM-Au₆, both the $\mathcal{w}$ value and the $\sigma$ parameter are lower than those in AzM-Au₄. In contrast, the hardness parameter decreases as the number of atoms in the gold cluster increases, with the $\eta$ value for AzM-Au₆ higher than that of AzM-Au₄. This observation suggests that the collective influence of the $\mathcal{w}$, ${\rm{ \mathsf{ σ} }}$, and $\eta$ parameters contributes to the heightened calculated Raman activity.

Table 4.  Some electronic properties for AzM and AzM-Au_n (n = 2, 4, and 6).

	ΔEbind (kcal/mol)	Ip (eV)	EA (eV)	χ (eV)
AzM	***	6.77	2.60	4.69
AzM-Au2	-21.67	5.70	2.93	4.31
AzM-Au4	-14.10	5.41	3.62	4.52
AzM-Au6	-15.10	5.87	2.97	4.42

---

 | Show Table

DownLoad: CSV

Table 5.  The reactivity parameters $\eta$, $\sigma$, and $\mathcal{w}$ (eV) for AzM and AzM-Au_n (n = 2, 4, and 6).

	$\eta$ (eV)	$\sigma$ (eV)-1	$\mathcal{w}$ (eV)
AzM	2.08426	0.47978	5.2785
AzM-Au2	1.38669	0.72114	6.7288
AzM-Au4	0.89675	1.11513	11.4061
AzM-Au6	1.44764	0.69077	6.77076

---

 | Show Table

DownLoad: CSV

3.1.5.   Electrostatic properties and molecular electrostatic potential (MEP) surfaces

Table 6 presents the reported electrostatic properties of the AzM-Au_n complexes. There is a consistent rise in the overall polarizability from AzM-Au₂ to AzM-Au₄. However, the value for AzM-Au₆ is lower compared to both AzM-Au₂ and AzM-Au₄. Consequently, we can deduce that the enhancement in Raman activity from AzM-Au₂ to AzM-Au₄ is directly associated with the augmentation in the total polarizability value. This is why AzM-Au₆ does not exhibit the anticipated rise in Raman activity.

Table 6.  Some electrostatic properties for AzM and AzM-Au_n (n = 2, 4, and 6).

	µ (Debye)	α0(a.u.)	Δα (a.u.)	|Δq| (electrons)
AzM	4.90	225.69
AzM-Au2	11.17	324.99	20.45	0.24
AzM-Au4	8.70	448.70	55.66	0.18
AzM-Au6	13.09	511.49	28.58	0.21

---

 | Show Table

DownLoad: CSV

The determination of MEP surfaces was performed using Gauss View 6 software. Illustrated in Figure 7 are the MEP surfaces corresponding to AzM and AzM-Au_n complexes. Examination of the figure discloses that within the pristine AzM molecule, the region exhibiting the highest electron density is in proximity to the sulfur atom, primarily attributed to the presence of lone pairs. Concurrently, there is a noticeable presence of delocalized electron density around the oxygen atom. Following the integration of silver clusters, a discernible overall reduction in electron density becomes apparent when compared to the unmodified AzM molecule. This reduction is attributed to the overlap of electron clouds with those originating from the gold clusters.

Figure 7.  The MESP surfaces of (a) AzM, (b) AzM-Au2, (c) AzM-Au4, and (d) AzM-Au6.

DownLoad: Full-Size Img PowerPoint

3.1.6.   Binding energies

Electronic energy was considered in the calculations of the binding energy (ΔE_bind). Although the configuration of AzM-Au₄ in Figure 5 demonstrates the anticipated rise in Raman activity, there is a departure from the trend observed in the reported binding Energies, as depicted in Table 4. As we know, a more stable configuration would manifest as an increasingly negative binding energy. However, the ΔE_bind for AzM-Au₄ does not adhere to this trend, suggesting that AzM-Au₄ is less stable when compared to AzM-Au₂. The heightened Raman intensity in AzM-Au₄ is associated with Global Electrophilicity Indices ($\mathcal{w}$), electronic propriety, and total polarizability (Δα). Conversely, the AzM-Au₆ complex is less stable than AzM-Au₂, as indicated by its higher value than AzM-Au₂. Additionally, AzM-Au₆ does not follow the expected increase in Raman intensity, aligning with trends observed in global electrophilicity indices ($\mathcal{w}$), electronic proprieties, and total polarizability (Δα), as detailed in the aforementioned tables.

3.1.7.   NBO and AIM analysis

We conducted an NBO Analysis on AzM-Au_n (n = 2, 4, and 6) complexes to delve deeper into the characteristics and durability of the interaction between the gold cluster and AzM. Table 6 shows the charge transfer values derived through the Natural Population Analysis (NPA) as an integral component of the NBO calculations. The Δq is characterized as the algebraic sum of the natural charges on the gold atoms. In Table 6, the magnitudes of the Δq for the AzM-Au_n complexes are presented, revealing a pattern akin to the trends observed in interaction energies and E² energies. This correlation underscores the connection between the stabilities of the complexes and the exchange of charges between the molecule and the gold cluster. Table 7 shows the E² second-order stabilization energies between sulfur and gold atoms for the AzM-Au_n (n = 2, 4, and 6) complexes. The trend of E² energies reflects the trend of the binding energies. The AzM-Au₂ complex is more stable due to a stronger overlap between the lone pairs of the two fragments. Compared to AzM-Au₄ and AzM-Au₆, both the E² energies and binding energies (without the minus sign) are lower, indicating their decreased stability compared to AzM-Au₂. In Table 8, we find some important parameters related to AIM analysis. The parameter ∇² ρ(r) is related to the energy associated with bond interactions. A positive value of ∇² ρ(r) indicates a decrease in charge density between bonds, while a negative value indicates an increase in charge density within the bond. Ionic interactions are associated with positive ∇² ρ(r) values, whereas covalent interactions are associated with negative values. The accumulation of charge at the (r) position results in H(r) taking on a negative value, indicating a stabilizing effect. When both the second derivative of the charge density concerning the position (∇² ρ(r)) and the Hamiltonian function (H(r)) are positive, the interaction is classified as an electrostatic bond. This bond can have varying strengths, ranging from weak to moderate hydrogen bonds. A polar covalent bond is identified by a positive ∇² ρ(r) value and a negative H(r) value, which means it is a combination of partial covalent and partial electrostatic characteristics. This type of bond is also referred to as a coordination bond or a strong hydrogen bond. If both atoms exhibit negativity, the bond is nonpolar covalent. Table 8 provides bond classifications according to AIM.

Table 7.  Results of the second-order perturbation theory analysis of the NBO basis of the optimized structures of AzM-Au_n complexes.

Complexes	Donor NBO (i)	Acceptor NBO (j)	E2 (kcal.mol−1)
AzM-Au2	ns2	n*(Au32)	24.33
AzM-Au4	ns2	n*(Au32)	13.66
AzM-Au6	ns2	n*(Au32)	14.60

---

 | Show Table

DownLoad: CSV

Table 8.  BCP data from the AIM Analysis of optimized AzM-Au_n complex structures.

Complexes	BCP	ρ(r)	∇2ρ(r)	H(r)	Interaction type
AzM-Au2	Au32−S	0.0707	0.1611	–0.1870	Polar covalent bond
	Au32····C-H	0.6122	0.0166	0.0006	Weak electrostatic bond
AzM-Au4	Au32−S	0.0620	0.1422	–0.0153	Polar covalent bond
	Au32····C-H	0.0066	0.0177	0.0007	Weak electrostatic bond
AzM-Au6	Au32−S	0.0620	0.1389	–0.0154	Polar covalent bond

---

 | Show Table

DownLoad: CSV

3.2.   Monte Carlo simulation

The adsorption energy in the context of Monte Carlo simulations refers to the amount of energy absorbed when a molecule is adsorbed onto a surface. In Monte Carlo simulations, which are a computational technique used to model complex systems, including molecular interactions, the adsorption energy is calculated as part of the simulation process. The adsorption energy provides insights into the stability and nature of the adsorption process. A more negative adsorption energy generally indicates a stronger adsorption interaction. It is an important quantity in understanding and predicting adsorption phenomena in various fields, including chemistry, materials science, and environmental science. The adsorption energies ${\varDelta E}_{ads}$ of AzM on diverse AuNP shapes (cubic box, right circular cone, cylinder, and sphere) and sizes (3, 5, and 7 nm) were thoroughly investigated using MC simulation and are summarized in Table 9. The lowest-energy structures of the adsorbed AuNP-AzM complexes are illustrated in Figure 8.

Table 9.  The adsorption energies (${\varDelta E}_{ads}, kcal.{mol}^{-1}$) values of AzM on various AuNP shapes.

${\varDelta E}_{ads}, kcal.{mol}^{-1}$
Shape	Sphere			Cubic			Cylinder			Right circular cone
Size	3nm	5nm	7nm	3nm	5nm	7nm	3nm	5nm	7nm	3nm	5nm	7nm
AzM	-84.14	-98.61	-102.78	-97.86	-98.11	-92.47	-97.78	-97.78	-97.71	-102.35	-105.32	-113.51

---

 | Show Table

DownLoad: CSV

Figure 8.  The lowest-energy structures of the AzM adsorption on 3nm, 5nm, and 7nm of Sphere, Cubic, Cylinder and Right circular cone AuNPs, as acquired through the MC simulation.

DownLoad: Full-Size Img PowerPoint

The impact of NP size on adsorption energies was evaluated, revealing that an increase in NP size has a negligible effect on the adsorption energies for the cylinder shape. Notably, for the cubic shape, the ${\varDelta E}_{ads}$ value at 7 nm is lower than those at 3 nm and 5 nm. Similarly, the ${\varDelta E}_{ads}$ values for the right circular cone and sphere shapes at 7 nm are greater than those at 3 nm and 5 nm. Interestingly, the ΔE_ads values for the 3 nm and 5 nm sizes (Right circular cone and cubic) are very close to each other. Based on our findings, it can be inferred that the right circular cone-shaped structure with a size of 7 nm exhibits a more pronounced negative adsorption energy, suggesting a stronger interaction in the adsorption process.

Table 10 provides a comparative analysis of various studies examining the interaction between molecules and nanomaterials using density functional theory (DFT) and Monte Carlo (MC) methods. Although there are no direct studies involving Azinphos-methyl, this comparison highlights similarities and differences in interaction behaviors across different molecules and nanomaterials. By analyzing how our findings align with existing research in this field, we offer a deeper understanding of the factors influencing molecular interactions with nanomaterials. This approach not only enhances the interpretation of our study but also situates it within the broader scope of molecular interaction research, offering a clearer understanding of how our work contributes to the field.

Table 10.  Comparative analysis of various studies examining the interaction between molecules and nanomaterials.

References	Complexes	Methods
[37]	C(AA)-Aun (n=2–10)	DFT-MC
[45]	C(GA)-Aun (n=2–10)	DFT-MC
[46]	Fipronil-Ag6	DFT
[47]	DMAB-Aun (n=3, 5, 19)	DFT
[36]	Tetrazole derivatives-Cu(111)	DFT-MC
[48]	pyrazolyl nucleosides-Cu(111)	DFT-MC
This work	AzM-Aun (n=2, 4, and 6)	DFT-MC

---

 | Show Table

DownLoad: CSV

4.   Conclusions

In conclusion, the theoretical investigation into the interaction of the Azinphos-methyl (AzM) molecule with gold nanoclusters and nanoparticles has yielded significant insights into the electronic properties, charge transfer mechanisms, and overall stability of these complexes. Utilizing a combination of density functional theory (DFT) calculations, natural bond orbital (NBO) analysis, and Monte Carlo (MC) simulations, the study offers a comprehensive understanding of how various factors such as the size and shape of gold nanoparticles influence the adsorption and interaction energies of Azinphos-methyl. Key findings indicate that the interaction between AzM and gold clusters substantially affects Raman activity and electronic properties. The increase in Raman intensity is closely associated with a rise in global electrophilicity indices ($\mathcal{w}$) and total polarizability (Δα), underscoring the importance of this interaction in potential applications. Computational analyses further reveal that as the number of gold atoms in the cluster increases, there is a notable enhancement in the stability of the AzM-Au_n complexes. This stability is correlated with increased charge transfer (Δq), binding energies (ΔE_bind), and E² energies, highlighting the critical role of cluster size in determining the strength of the interaction. However, among the complexes studied—AzM-Au₂, AzM-Au₄, and AzM-Au₆—the AzM-Au₂ complex stands out as the most stable. This finding suggests that smaller gold clusters may provide optimal conditions for stable complex formation with Azinphos-methyl. The Monte Carlo simulations offer additional insights into the adsorption process, indicating that a right circular cone-shaped gold nanoparticle, particularly with a size of 7 nm, exhibits the most substantial negative adsorption energy. This suggests a stronger interaction during the adsorption process, causing this geometry and size to be favorable for the adsorption of AzM. In summary, this study not only advances our understanding of the interaction between Azinphos-methyl and gold nanostructures but also provides a foundation for future research and potential applications in areas such as sensor development, environmental monitoring, and nanotechnology-based pesticide detection and removal.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors extend their appreciation to the Researchers Supporting Project (RSPD2024R808) of King Saud University, Riyadh, Saudi Arabia.

Conflict of interest

The authors declare no conflict of interest.