Adsorption of non-steroidal anti-inflammatory drugs (NSAIDs) on nanographene surface: Density functional theory study

1 Introduction

The existence of pharmaceutical pollution resulting from human use in the environment has been a major topic for over 25 years now (Kümmerer, 2009). Hundreds of tons of pharmaceuticals are being produced, dispensed, and consumed annually. After being discarded by patients or health care professionals, many drugs and their metabolites make their way into the environment (Sherer, 2006). Therefore, they have become universal soil and water pollutants that may have permanent detrimental effects on aquatic organisms and human health. Pharmaceuticals or active pharmaceutical ingredients (APIs) are complex molecules with different functionalities and physio-chemical and biological properties. They are developed and used because of their specific biological activity. Most of them are polar compounds (Kümmerer, 2009).

There has been a rapid growth in research that has focused on the environmental removal of non-steroidal anti-inflammatory drugs (NSAIDs). This is due to their wide usage for the treatment of numerous diseases such as Musculo-skeletal pain (D’Angelo et al., 2012; Alturki et al., 2018), arthritis (Sultana and Rasool, 2015), and other pains (Asano et al., 2019; Planas et al., 2003), and a patient can take them without a medical prescription (Mlunguza et al., 2019). This leads to high concentrations (in μg/L) of these drugs in aquatic environments worldwide. Poor control of NSAIDs in the environment can result in their uptake by crops through farming activities.

NSAIDs are weak organic acids having a carboxylate unit (Mlunguza et al., 2019) with a pK_a of 3-to-5. There are many NSAIDs drugs. The most used is Aspirin and many others. Necessary NSAIDs drugs are the targeted investigated drugs in this study and will be defined later in this section.

The difficulty of removal efficiency of NSAIDs drugs from household water and wastewater treatment plants (WWTPs) can be attributed to their: (1) higher solubility in water and (2) polar nature. One of the most widely used techniques for effective removal of NSAIDs is through adsorption. There are several adsorbents in plenty of research articles that have been proposed for the effective experimental removal of NSAIDs from the aquatic environment. Activated carbon (AC) in its two forms: (1) powdered activated carbon (PAC) and (2) granulated activated carbon (GAC) has been widely explored for the adsorption of NSAIDs. Luo et al. (2014) Other used adsorbents are the molecular imprinting polymer (MIP) (Samah et al., 2018), iron nano-adsorbent (Ali et al., 2016), metal-organic frameworks (MOFs) (Hasan et al., 2013).

Among the most important adsorbent system is graphene. Graphene can be modified into different forms. Therefore, it is used widely as an adsorbent of water contaminants. Additionally, graphene has a high surface area, which means it can extract more compounds and is not easily saturated by pollutants. This means that graphene-based adsorbents can adsorb plenty of pollutants at once (Liu et al., 2016). Moreover, nanoparticles of graphene, i.e., nanographene (NG) platelets, perform better as adsorbents because of their: (1) efficient working principle, (2) excellent optical and electrical properties, and (3) good chemical/physical stability (Hiew et al., 2018; Al-Khateeb et al., 2017; Zambianchi et al., 2017). NGs are known for their high surface areas, which result in high adsorption capabilities.

Al-khateeb and coworkers investigated experimentally (Al-Khateeb et al., 2017) the elimination of some NSAIDs, namely (3-(3-hydroxybut-3-en-2-yl)phenyl)(phenyl)methanone (Ketoprofen), (S)-2-(6-methoxynaphthalen-2 yl)propanoic acid (Naproxen), 3-(2-((2,6-dichlorophenyl)amino)phenyl)prop-1-en-2-olate (Diclofenac sodium salt), and 2-(4-isobutylphenyl)propanoic acid (Ibuprofen), abbreviated as KETO, NAP, DIC, and IBU, respectively, from an aqueous model and real solutions using NG surfaces, see Fig. 2 (a). They found that the % removal of the drugs from water by the NG surface is efficient and increased as the NG's concentration increased, see Fig. 1. In Fig. 1, data were collected from ref. Al-Khateeb et al. (2017) and represented graphically by the present author. Their results also showed that the values of Gibbs free energy, entropy, and enthalpy were negative, positive, and positive, respectively. These observations are evident of the spontaneity of the elimination process.

Close

Fig. 1

% removal of the four investigated NSAIDs from aqueous solutions by NG plates. Data were obtained from ref. (Al-Khateeb et al., 2017).

Export to PPT

Close

Fig. 2

(a) Chemical structures along with names and abbreviations, (b) optimized geometries of the four investigated NSAIDs, and (c) optimized geometry of NG surface.

Export to PPT

Quantum mechanical procedures especially DFT calculations, have proven to be an efficient method to study the interaction between nanomaterials and drugs (Wazzan et al., 2019). In comparison to the amount of experimental work carried in this regard, few theoretical studies were conducted. Kaczmarek-Kędziera studied NSAID diclofenac interactions with models of pristine chitosan and its modified chains (Singh et al., 2014) theoretically. Kaczmarek-Kędziera concluded from the interaction energy values of chitosan-diclofenac complexes and two chitosan chains that penetrate the drug molecule the chitosan surface is brutal, and the presence of the hydrogen bond donors and acceptors both in biopolymer and in diclofenac may cause the process to occur. Also, another finding in that study is regarding the modification of chitosan surface with long-distanced amino groups resulted in blocking the interaction between the two monomers (diclofenac and chitosan surface), keeping the interaction energy in the range from 20 to 40 kcal/mol. Singh and coworkers (Singh et al., 2014), investigated using quantum chemical calculations five NSAIDs, i.e., the acetic acid derivative of COX inhibitors. The Five NSAIDs are diclofenac, etodolac, indomethacin, ketorolac, and nabumetone. The study used the Mulliken atomic charges on each atom for each drug, the HOMO and LUMO energies and distributions to assist and interpret the results regarding the molecular interactions and the intramolecular interactions molecular packing of the five NSAIDs.

As an essential and complementary matter, theoretical study of the four drugs (KETO, NAP, DIC, and IBU) was carried out in their free and adsorbed forms on the NG surface for the first time. The electronic, electrostatic properties and some important quantum chemical parameters of the free drug molecules were calculated and compared. More importantly, the binding energies and thermodynamic parameters of the adsorbed systems of the drug over NG (drug@NG) were calculated and analyzed too. The calculated data were compared with the available experimental ones (Al-Khateeb et al., 2017).

2 DFT details

Gaussian 09 package (Frisch, 2016) was used to carry out the calculations. Visual inspections were performed using GaussView program (version 5.0.8) (Dennington et al., 2009) and Chemcraft program version 1.8 (build 489) (Zhurko and Zhurko, 2009). DFT/B3LYP (Wazzan et al., 2019; Shayan and Nowroozi, 2018; Hazrati et al., 2017) calculations were used in order to investigate the interaction between the nanographene (NG) and the four drugs IBU, KETO, DIC, and NAP molecules. DFT calculations with the reliable exchange-correlationB3LYP functional were performed. The B3LYP is a consistent functional which has been usually used in the study of similar nanostructures (Shayan and Nowroozi, 2018; Hazrati et al., 2017; Hosseinian et al., 2017; Baei et al., 2017; Bezi Javan et al., 2016; Baei et al., 2014). B3LYP functional is associated with the dispersion correction term “D” of Grimme (Becke, 1993; Lee et al., 1988; Hazrati et al., 2017) using Gaussian semiempirical dispersion = gd3 keyword. Indeed, most density functionals were incapable of describing the van der Waals interactions properly due to dynamical correlations between fluctuating charge distributions. A practical method to solve this issue was achieved by the DFT-D approach, which added a semiempirical dispersion potential to the conventional Kohn–Sham DFT energy (Duverger et al., 2014). The standard split-valance triple zeta polarized basis set (6-311G(d)) was used in all calculations (Ditchfield et al., 1971). The nanographene as an allotrope of carbon consisting of a single layer of atoms arranged in a two-dimensional honeycomb lattice surface. The modulization of NG surface in this study was performed using 28 hexagonal lattices consisting of 106 atoms (80 Carbon atoms saturated with 26 terminal Hydrogen atoms), Fig. 2 (c) shows the optimized geometry of NG surface.

For the adsorbed systems (drug@NG), the basis set superposition error (BSSE) was applied using the counterpoise (CP) method (Boys and Bernardi, 1970). BSSE typically originated due to the incompleteness description of the basis set in the calculations of week intermolecular interactions. (Abdolahi et al., 2018) The BSSE corrected adsorption energies ( $E_{a\mathrm{d}\mathrm{s}}^{\mathrm{B}\mathrm{S}\mathrm{S}\mathrm{E}}$ ) of all optimized systems were calculated using the following equation:

3 DFT results

3.1

3.1 Free drugs

The geometries of the four investigated NSAIDs were optimized at the B3LYP/6-311G(d) level. Each optimized geometry was confirmed to be a minimum point in the PES using the frequency calculations since no imaginary frequency (negative frequency) was obtained. The optimized geometries are shown in Fig. 2 (b).

Among the four optimized geometries, NAP and IBU show the most planar geometries. The two attached rings in NAP and the one ring in IBU allow a large part of the molecules' skeleton to lei in the same plane. On the other hand, the KETO and DIC geometries are significantly deviated from planarity, and the reason is the electrostatic repulsion between the two adjacent aromatic rings.

From the geometry optimization calculations, it is straightforward to obtain the frontier molecular orbitals (FMOs) energies and distributions. The analysis of FMOs is necessary in order to get valuable information about the molecular reactivity (Fukui et al., 1952). The energies and distribution illustration of the HOMO and LUMO of the four investigated drugs are shown in Fig. 2. The HOMO energy levels ( $E_{\mathit{HOMO}})$ is a measurement of how much the molecule can donate electron/s. As the $E_{\mathit{HOMO}}$ values are less negative (destabilized) as the molecule can donate electron/s more easily. The $E_{\mathit{HOMO}}$ values of the four investigated molecules increase in the following order: KETO < IBU < NAP < DIC. On the other hand, the LUMO energies ( $E_{\mathit{LUMO}})$ is a measurement of the ability of the molecule to receive an electron/s. It became destabilized (less negative) as the molecule became more able to accept electron/s. The order $E_{\mathit{LUMO}}$ values of the four investigated molecules increase in the order: KETO < NAP < DIC < IBU. Therefore, the KETO drug shows the least ability to donate and the largest ability to accept electron/s. While the IBU drug shows a moderate ability to donate and the least ability to accept electron/s.

According to Fig. 3, the HOMO and LUMO energies of the four NSAIDs are significantly different than those belonging to the NG plates. Therefore, it is expected for there to be a significant charge transfer between NSAIDs and NG. This agreed with the amount of charge transfer from the net natural atomic charges analysis, see Section 3.2.1. This type of interaction is ideal for the removal efficiency because the NSAIDs are expected to hold tightly to the NG surface without the possibility of easy release (Dastani et al., 2020).

Close

Fig. 3

HOMO and LUMO energy levels, and the numbers on the arrows are their energy gaps of the four investigated NSAIDs and NG, and the 3D distributions of HOMO and LUMO orbitals of the four investigated NSAIDs.

Export to PPT

The molecular reactivity and kinetic stability can be estimated by evaluating the energy gap. The energy gap ( ${\mathrm{\Delta }}_{H-L})$ is the difference in energies of the HOMO and LUMO orbitals. Generally speaking, a lower ${\mathrm{\Delta }}_{H-L}$ value indicates a higher molecular reactivity and lower kinetic stability and vice versa. The order of decreasing the ${\mathrm{\Delta }}_{H-L}$ values of the four investigated drugs is: IBU > KETO > NAP > DIC.

Fig. 3 shows the special distribution of the HOMO and LUMO orbitals over the skeleton of the four investigated drugs. The HOMO and LUMO orbitals are distributed over the whole skeleton of NAP, KETO, and IBU, especially in the case of NAP and IBU. However, the LUMO orbital is not delocalized over the carboxylic terminal in the KETO. The HOMO orbitals of DIC drug is distributed among most parts of the molecule. The Na atom is attached to the carboxylic terminal, while the LUMO orbital is mainly distributed on this atom (Na) and part of the carboxylic terminal. Fig. 4 represented the distribution of the total density mapped with the electrostatic potential abbreviated as MESP. The color ranges from red to blue indicating region where the electron density gradually decreases.

Close

Fig. 4

Electrostatic properties of the four investigated NSAIDs represented by 3D-MESP surfaces and 2D-contour maps.

Export to PPT

The quantum chemical parameters (QCPs) were calculated at the B3LYP/6-311G(d) level and listed in Table 1. The calculated QCPs included the ionization potential ( $\mathit{IP}$ ), electron affinity ( $\mathit{EA}$ ), dipole moment ( $\mu$ ), electronegativity ( $\chi$ ), hardness ( $\eta$ ), softness ( $\sigma$ ), and electrophilicity index ( $\omega$ ) (Parr, 1999; Pearson, 1986). The equation beside each parameter in Table 1 indicates the method used to evaluate such parameters.

Table 1 QCPs of the four investigated NSAIDs.

QCPs	KETO	DIC	NAP	IBU
Ionization potential $(IP=-E_{\mathit{HOMO}})$	7.17	5.47	5.70	6.60
Electron affinity $(EA=E_{\mathit{LUMO}})$	2.29	1.59	1.26	0.53
Dipole moment $(\mu )$	2.92	5.59	1.85	1.37
Electrophilicity index $\left(\omega =\frac{{\eta }^2}{2\eta }\right)$ (Parr, 1999)	47.45	219.39	21.05	8.44
Global hardness $\left(\eta =\frac{\mathit{IP}-EA}{2}\right)$ (Parr and Pearson, 1983)	2.44	1.94	2.22	3.03
Global softness $\left(\sigma =\frac{1}{\eta }\right)$ (Parr and Pearson, 1983)	0.41	0.52	0.45	0.33
Electronegativity $\left(\chi =\frac{\mathit{IP}+EA}{2}\right)$ (Parr and Pearson, 1983)	4.73	3.53	3.48	3.56

Regarding studying the role each molecule can play, i.e., which molecule is an electron-donor, an electron-acceptor, or an electron-donor and electron-acceptor simultaneously, it is very convenient to analyze the $\mathit{IP}$ and $\mathit{EA}$ as these two quantities measure these two abilities. The $\mathit{IP}$ is the minimum amount of energy required to remove the valence electron of an isolated neutral gaseous atom or molecule. The order of increasing the $\mathit{IP}$ values and decreasing the ability to act as electron-donor is: DIC < NAP < IBU < KETO. On the other hand, the $\mathit{EA}$ is the amount of energy released when an electron is added to a neutral atom or molecule in the gaseous state to form a negative ion. The order of increasing the $\mathit{EA}$ values and increasing the ability to act as electron-acceptor is: IBU < NAP < DIC < KETO.

The two important global parameters that can be used to discuss the reactivity of a drug toward adsorption on the NGs surface is the chemical hardness ( $\eta$ ) and softness ( $\sigma$ ) parameters. The values of $\eta$ and $\sigma$ are frequently associated with the Pearson’s hard and soft acids and bases (HSAB) (Pearson, 1986). From the obtained results, the hard molecule has a large ${\mathrm{\Delta }}_{H-L}$ and therefore is less reactive, whereas the soft molecule has a smaller ${\mathrm{\Delta }}_{H-L}$ and is therefore more reactive. The order of decreasing the $\eta$ values is the same as that of ${\mathrm{\Delta }}_{H-L}$ values, i.e., IBU > KETO > NAP > DIC.

Higher electronegativity values show that the investigated drug shows high reactivity and has the propensity to share more electrons involved in the chemical bond. The values of $\chi$ of the four investigated drugs are comparable. The difference between the maximum and minimum value is only 1.25 eV. However, the order of decreasing the $\chi$ values is: KETO > IBU $\ge$ DIC > NAP. Therefore, the KETO drug shows the largest $\chi$ value and is expected to be the most reactive molecule as an electron-acceptor. Similarly, the electrophilicity index ( $\omega$ ) is a measure of the ability of a molecule to receive electrons, and therefore, for a molecule, it measures the electrophilicity ability. The $\omega$ values of the four investigated drugs decrease in the following order: DIC > KETO > NAP > IBU. Therefore, the DIC molecule with a significant larger $\omega$ value, because of the sodium acetate terminal, is expected to be the most reactive molecule as an electron-acceptor.

As reported, the adsorption process is influenced by hydrogen bonding (H-Bonding), van der Waals interactions, and π-π interactions (Al-Khateeb et al., 2017; Baccar et al., 2012), since the presence of the π electrons on both organic compounds and the NG surface can lead to the π-π layering of the drugs on the NGs (Al-Khateeb et al., 2017). On the other hand, it is expected that the dipole-dipole interactions between the NGs surface and the adsorbed drug will be enhanced by the larger dipole moment ( $\mu$ ) value, which result in better adsorption. The $\mu$ values of the investigated drug molecules, Table 1, increase in the order: IBU < NAP < KETO < DIC. The dipole moments of the KETO and DIC drugs are large compared to that of water (1.8546 Debye). This result reveals that the dipole-dipole interactions between these two drugs and the NGs are significant.

3.2

3.2 Drug@NG adsorbed systems

To achieve the π-π stacking configuration between the drug and the adsorbent, each drug molecule is initially placed horizontally above the NG surface, and the geometry optimization process of the adsorbed system (drug@NG) was carried out free of any structural constraints. The horizontal orientation was chosen to achieve the π-π stacking between the π-electrons of the drug molecule and the NGs. After the completion of the optimization process, the drug molecule oriented itself in a vertical position above the NG surface, and the exception is DIC@NG system, where the stacked orientation is almost kept. The optimized systems of drug@NG are presented in Fig. 5. However, this is expected for the NAP drug due to its mostly planetary structure, followed by IBU adsorbed in the π-π stacking configuration. IBU@NG is almost in the π-π stacking configuration, while the NAP@NG shows the clearer vertical (T-shaped) orientation. In this case we could attribute the vertical positioning of NAP drug during the adsorption process to the stability formed due to H-bonding formation. Thus, for this particular drug, the H-bonding formation's stability is overcome due to the π-π stacking configuration (Hasanzade and Raissi, 2019).

Close

Fig. 5

Optimized systems of NSAID@NG adsorbed systems.

Export to PPT

3.2.1

3.2.1 Natural population analysis

Fig. 6 shows the plotting of the net natural atomic charges ( $\sum$ NPA charges/e) of the atoms belonging to the two components of the adsorption systems (KETO, DIC, NAP, IBU component, NG component) as obtained from the B3LYP-D/6-311G(d) level. The calculation of NPA charges only involves calculating the natural atomic orbitals and summation over all-natural atomic orbital occupancies (NAOs) of a given atom to obtain the natural charges for each of the atoms.

Close

Fig. 6

Plotting of natural atomic charges of the two components of drug@NG.

Export to PPT

Our results showed that the $\sum$ NPA on the isolated components are zero. After the adsorption process, $\sum$ NPA values become either a negative or positive value. When $\sum$ NPA values are negative for one component, this indicates that the electrons are transferred to it from the other component during the adsorption process and vice versa.

It is worth mentioning that the negative and positive values of $\sum$ NPA are quite small, indicating that the electrostatic interaction between the two components of the adsorbed systems is weak. However, in two systems, i.e., DIC@NG and IBU@NG, the $\sum$ NPA values of the drug component are negative, and they are $-0.0120$ and $-0.0060$ e, respectively. On the other hand, the NG components have positive and equal values of $\sum$ NPA. Therefore, in these two systems, the drug acts as an electron-acceptor, and the NG acts as an electron-donor. For the other two adsorbed systems, i.e., KETO@NG and NAP@NG, the reverse is true. The KETO and NAP drug components carry the partially positive charges, 0.0010 and 0.0002 e, respectively, and the NG components carry the partially equal and negative charges. The KETO and NAP are expected to act as electron-donors, and the NG acts as electron-acceptor. An interesting notice is that the amount of charge transfers in the DIC@NG is significantly larger than those in the other three systems; this will undoubtedly expect to enhance the resultant adsorbed system's stability, as evident experimentally (Al-Khateeb et al., 2017).

3.2.2

3.2.2 Adsorption energies and thermodynamic parameters

Fig. 7 shows the BSSE and the BSSE adsorption energies ( ${\mathrm{\Delta }E}_{\mathit{ads}}^{\mathit{BSSE}}$ ) values for the four adsorbed systems, IBU@NG, DIC@NG, KETO@NG, and NAP@NG, as calculated at the B3LYP-D/6-311G(d) level. The ${\mathrm{\Delta }E}_{\mathit{ads}}^{\mathit{BSSE}}$ values are less than 3.89 kcal/mol. The BSSE values are 2.96, 3.89, 2.85, and 2.27 kcal/mol for complexes IBU@NG, DIC@NG, KETO@NG, and NAP@NG, respectively. The relatively larger values of BSSE indicates the importance of including this correction in the calculations of the binding energies (Al-Qurashi and Wazzan, 2017).

Close

Fig. 7

The BSSE corrected adsorption energies ( $\mathrm{\Delta }{E}_{\mathit{ads}}^{\mathit{BSSE}}$ ), and (−) BSSE energy ( $-E^{\mathit{BSSE}}$ )(in kcal/mol) of the adsorbed systems (drug@NG).

Export to PPT

As the binding energy became more negative, this indicated a tide interaction between the drug and nanographene surface. The binding energies decrease (become more negative) in the following order: NAP@NG (−109.23 kcal/mol) > IBU@NG (−115.74 kcal/mol) > KETO@NG (−118.42 kcal/mol) > DIC@NG (−127.56 kcal/mol). The values of the adsorption energies are quite large and may signify a chemisorption process. The theoretical finding regarding the larger binding energy of the DIC@NG system compared to other adsorbed systems matches the experimental findings (Al-Khateeb et al., 2017). In the experimental testing, it was found that %removal of DIC drug from water by NG is the maximum. However, the experimental % removal with 10 mg of NG of the four drugs from aqueous solutions was: IBU (62.9%) < KETO (83.3%) < NAP (91.2%) < DIC (95.4%), not very consistent with the trend observed for the binding energies.

The thermochemistry of the adsorption process has been obtained. The obtained thermochemical parameters included the sum of electronic energy corrected with zero-point energy changes $({\mathrm{\Delta }E}_{\mathit{ele}}^{\mathit{ZPE}})$ , enthalpy changes $({\mathrm{\Delta }H}_{\mathit{ads}}$ ), Gibbs free energy changes ( ${\mathrm{\Delta }G}_{\mathit{ads}}$ ) and entropy changes ( ${\mathrm{\Delta }S}_{\mathit{ads}}$ ) are collected in Table 2.

Table 2 Thermochemical parameters for the adsorption systems of the four drugs on NG surface at the 298.15 K calculated at the B3LYP-D/6-311G(d) level of theory. Note: in parenthesis are experimental values from ref. (Al-Khateeb et al., 2017).

Adsorption system	${\mathrm{\Delta }E}_{\mathit{ele}}^{\mathit{ZPE}}$ /kcal mol−1	${\mathrm{\Delta }H}_{\mathit{ads}}$ /kcal mol−1	${\mathrm{\Delta }G}_{\mathit{ads}}$ /kcal mol−1	${\mathrm{\Delta }S}_{\mathit{ads}}$ /cal mol−1.K−1	$\mathit{ZPE}$ /kcal mol−1
KETO@NG	−116.67	−116.01 (+9.4)	−105.02 (−0.81)	−36.84 (+34.60)	653.62
DIC@NG	−126.31	−121.13 (+11.2)	−113.91 (−1.78)	−24.21 (+43.97)	620.87
NAP@NG	−107.96	−107.20 (+10.2)	−98.12 (−1.14)	−30.47 (+38.17)	646.28
IBU@NG	−115.24	−113.51 (+6.57)	−102.88 (−0.29)	−35.67 (+23.21)	667.72

From a thermodynamic perspective, the adsorption process is favorable, since the values of ${\mathrm{\Delta }H}_{\mathit{ads}}$ , ${\mathrm{\Delta }G}_{\mathit{ads}}$ , and ${\mathrm{\Delta }S}_{\mathit{ads}}$ are all negative. The values of ${\mathrm{\Delta }H}_{\mathit{ads}}$ become more negative following the order: NAP@NG (-107.20 kcal/mol) > IBU@NG (-113.51 kcal/mol) > KETO@NG (-116.01 kcal/mol) > DIC@NG (-121.13 kcal/mol). The values of ${\mathrm{\Delta }G}_{\mathit{ads}}$ are becoming more negative following the same order observed for the ${\mathrm{\Delta }H}_{\mathit{ads}}$ values, and they are: −98.12, −102.88, −113.91, 105.02, and −113.91 kcal/mol for NAP@NG, IBU@NG, KETO@NG, and DIC@NG, respectively. In contrast, the values of ${\mathrm{\Delta }S}_{\mathit{ads}}$ become more negative following a different order, i.e., DIC@NG (-24.21 cal/mol.K) > NAP@NG (-30.47 cal/mol.K) > IBU@NG (-35.67 cal/mol.K) > KETO@NG (-36.84 cal/mol.K). Therefore, for the four investigated drugs, the order of increasing the ${\mathrm{\Delta }E}_{\mathit{ads}}^{\mathit{BSSE}}$ values (become more negative) is following the same order as the values of ${\mathrm{\Delta }H}_{\mathit{ads}}$ and ${\mathrm{\Delta }G}_{\mathit{ads}}$ . This result indicates two points: (1) as the ${\mathrm{\Delta }H}_{\mathit{ads}}$ and ${\mathrm{\Delta }G}_{\mathit{ads}}$ values become more negative as the interaction between the drug, and the NG surface becomes more spontaneous, and (2) the adsorption of a drug is enhanced as the drug can be adsorbed more spontaneously on the NG surface. The order in which the ${\mathrm{\Delta }S}_{\mathit{ads}}$ values become more negative is not consistent with the order of increasing the ${\mathrm{\Delta }E}_{\mathit{ads}}^{\mathit{BSSE}}$ values. However, DIC@NG and NAP@NG adsorbed systems with the maximum and minimum ${\mathrm{\Delta }E}_{\mathit{ads}}^{\mathit{BSSE}}$ values, respectively, have the largest and smallest ${\mathrm{\Delta }S}_{\mathit{ads}}$ values, respectively. It is maybe suitable to assume that the adsorption of the drug is stronger as the randomness (disorder) of the resultant system is decreased and consequently the ${\mathrm{\Delta }S}_{\mathit{ads}}$ values become less negative. This is confirmed by the observation that the ${\mathrm{\Delta }E}_{\mathit{ele}}^{\mathit{ZPE}}$ values are more negative compared to the ${\mathrm{\Delta }G}_{\mathit{ads}}$ values which shows that due to the adsorption process, the entropy is increased due to the decrease of disorder of the system (Wazzan et al., 2019; Al-Qurashi and Wazzan, 2017). However, in comparison to the experimental data in ref. (Al-Khateeb et al., 2017), the values of ${\mathrm{\Delta }H}_{\mathit{ads}}$ , ${\mathrm{\Delta }G}_{\mathit{ads}}$ , and ${\mathrm{\Delta }S}_{\mathit{ads}}$ are positive, negative, positive, respectively. Therefore, there is an agreement between the calculated and experimental negative values of ${\mathrm{\Delta }G}_{\mathit{ads}}$ . Experimentally and theoretically, the DIC drug adsorbed at the NG surface shows the most negative ${\mathrm{\Delta }G}_{\mathit{ads}}$ value, indicating that the adsorption is the most spontaneous one (see Table 2).

4 Conclusions

The adsorption of the four NSAIDs on the NG surface as removal was investigated at the B3LYP-D/6-311G(d) level of theory.

• The three drugs reoriented themselves during the optimization process in a vertical position except for the DIC drug.

• The adsorption energies of all investigated systems are negative; negative binding energies indicate the stability of the formed adsorbed system relative to its constituent monomers.

• Due to the relatively larger $E_{\mathit{ad}}^{\mathit{BSSE}}$ for the four adsorbed systems, theoretically, it is expected that the adsorption of these drugs on the NG surface is through a chemisorption process, and this gives an advantage of NG as a removal for the four drugs.

• The adsorption of the four drugs is arranged according to the increase of $E_{\mathit{ad}}^{\mathit{BSSE}}$ values and consequently the stability of the formed adsorbed system as: NAP@NG < IBU@NG < KETO@NG < DIC@NG.

• However, the order of the calculated stabilities of the formed complexes is not agreed with the order of the % removal reported experimentally, indicating the impact of other factors upon the adsorption process such the pH, solubility etc.

• Among the four systems, the DIC@NG system showed the largest binding energy (most stable system), totally agreed with the experimental finding.

• Computed negative free energy changes matched the negative experimental values but did not agree with the reported trend or values. However, the DIC@NG, the most stable adsorbed system, and the system with the highest % removal value shows the most negative ${\mathrm{\Delta }G}_{\mathit{ads}}$ .

• The spontaneity of the adsorption process was advent from the negative values of the thermodynamic parameters.