Anticorrosive Effects of Some Thiophene Derivatives Against the Corrosion of Iron: A Computational Study

It is known that iron is one of the most widely used metals in industrial production. In this work, the inhibition performances of three thiophene derivatives on the corrosion of iron were investigated in the light of several theoretical approaches. In the section including DFT calculations, several global reactivity descriptors such as EHOMO, ELUMO, ionization energy (I), electron affinity (A), HOMO-LUMO energy gap (ΔE), chemical hardness (η), softness (σ), as well as local reactivity descriptors like Fukui indices, local softness, and local electrophilicity were considered and discussed. The adsorption behaviors of considered thiophene derivatives on Fe(110) surface were investigated using molecular dynamics simulation approach. To determine the most active corrosion inhibitor among studied thiophene derivatives, we used the principle component analysis (PCA) and agglomerative hierarchical cluster analysis (AHCA). Accordingly, all data obtained using various theoretical calculation techniques are consistent with experiments.


INTRODUCTION
One of the serious problems in industrial sector is the corrosion of metals or alloys, which causes great casualties and enormous property loss (Frankel et al., 2015;Li et al., 2015). The most environmentally-benign and cost-effective approach to prevent metals against corrosion in acid solutions is using of inhibitors (Raja et al., 2016). Organic molecules containing O, N, and/or S atoms are the most widely used and are considered to be effective corrosion inhibitors (Xhanari et al., 2017). It is generally assumed that they can be adsorbed at metal surface through some active groups like heteroatoms, triple bonds or aromatic rings (Kovačević and Kokalj, 2013;Ko and Sharma, 2017). Previous studies showed that most organic inhibitors decrease corrosion rate by adsorption on the substrate surface and the inhibition performance follows the sequence O < N < S < P (Loto et al., 2012). Bockris and Swinkels suggested that S and/or N atoms could easily adsorb on metal surface by replacing water molecules (Bockris and Swinkels, 1964). The interactions between the organic corrosion inhibitors and the metal substrates are generally divided into two types: physical adsorption and chemical adsorption. On the basis of the theory hard and soft acids and bases (HSAB) introduced by Pearson (Pearson and Songstad, 1967), molecules including sulfur atom in their molecular structures can be regarded as soft bases, which give easily electron to metal surfaces because they cannot resist to electron cloud polarization or deformation. Thiophene is a sulfur containing heterocyclic compound with the molecular formula C 4 H 4 S. The electron pairs on sulfur atom are delocalized in π-conjugated systems. Thiophene and its derivatives can be obtained from petroleum or coal and in general are wellknown because of their therapeutic applications in medicinal chemistry. Many theoretical and experimental studies including the analyzing of inhibition efficiencies of thiophene and its derivatives are available in the literature (Benabdellah et al., 2011;Gece, 2013;Yadav et al., 2014). Especially, Fouda et al. (2014) synthesized three thiophene derivatives and investigated their anticorrosive effects against the corrosion of carbon steel using experimental methods such as weight loss, Tafel test, electrochemical frequency modulation, and electrochemistry impedance test. The molecular structures of mentioned molecules are given in Figure 1.
It is well-known that reactivity in chemistry is a backbone because it is well associated with reaction mechanisms. Therefore it allows us to understand chemical reactions, behavior of substance and improving synthesis procedures to obtain new materials such as corrosion inhibitors or drugs. Quantum chemistry calculations have proven to be very effective in evaluating the corrosion inhibition efficiency (Khalil, 2003;Obot et al., 2015;Taylor, 2015;Lgaz et al., 2018). They have been widely employed to interpret the experimental phenomena. Especially in recent years, conceptual density functional theory (CDFT) has been developed and applied to analyze the molecular activity of inhibitors (Geerlings et al., 2003;Liu, 2009). Furthermore, molecular dynamics simulation has become another effective way to explore the structure and bonding characteristics at inhibitor/metal interface (Khaled, 2008;Oguzie et al., 2013;Wang et al., 2016).

Gaussian Calculations
Gaussian 09 package program was used to the study of the isolated compounds using density functional theory (DFT), with the B3LYP functional (Wiberg, 2004) under the SDD, 6-31G and 6-31++G basis sets. For comparison, the calculations based on Hartree-Fock (HF) theory were also performed. Gauss View 5.0.8 program was employed to prepare the correlative calculation parameters for the compounds under probe.
As we all know, electrochemical corrosion generally happens in liquid environment. One of the most popular approaches to research solvent effect is to consider hydrogen bonded clusters of solvent molecules surrounding the solute molecules. Thus selfconsistent reaction field (SCRF) theory, with Tomas's polarized continuum model (PCM) (Andzelm et al., 1995) was used to describe the solvent effect of water. This method describes the solvent as a structureless continuum with uniform dielectric permittivity, in which a molecular-shaped empty cavity is dug to host the solute (Scalmani et al., 2006;Aouniti et al., 2013). The reliability of PCM method to explore the solvent effect in the field of corrosion inhibitors has been validated by many researchers (Wazzan et al., 2016;Guo et al., 2017;Yang et al., 2017).

Global Reactivity Descriptors
Chemical reactivity can be simply defined as the tendency of a chemical matter to undergo chemical reaction with another chemical matter. It is well-known that the understanding of the nature of chemical interactions and the prediction of chemical reactivity of atoms, ions or molecules are some of the challenging issues in chemistry. In the CDFT, quantum chemical descriptors like electronegativity (χ), chemical hardness (η), and chemical potential (µ) are usually considered. µ and η are defined as the first derivative of the electronic energy and chemical potential with respect to the electron number (N) at constant external potential, v (r) , respectively (Liu, 2009;Frau and Glossman-Mitnik, 2017).
Within the framework of finite differences approximation, the following expressions based on first vertical ionization energy and electron affinity values of chemical species are given (Madkour and Elshamy, 2016).
Softness that is known as a measure of the polarizability is mathematically defined as the multiplicative inverse of chemical hardness: Ionization energies and electron affinities of molecules can be predicted via Koopman's Theorem (Bellafont et al., 2015), i.e., ionization energy and electron affinity values of a chemical species correspond to negative values of its HOMO and LUMO orbital energies, respectively. So we can write the following equations for the calculation of quantum chemical parameters like hardness, electronegativity and chemical potential.
Besides, global electrophilicity index (ω) introduced by Parr (Parr et al., 1999), nucleophilicity (ε), which is physically the inverse of the electrophilicity, the fraction of electrons ( N) transferred from (to) the inhibitor molecule to (from) the metal surface, the energy of back donation E b-d , electronic charge accepting capability and the initial molecule-metal interaction energy ψ, were calculated in terms of global hardness (η) and electronegativity (χ) as given in .  Herein, χ Fe , χ inh , η Fe , and η inh represent the absolute electronegativity and hardness of iron and inhibitor molecule, correspondingly. In order to obtain the N values, we adopted a theoretical value of χ Fe = 7.0 eV and η Fe = 0 by assuming that for a metallic bulk I = A since they are softer than the neutral metallic atoms (Zarrouk et al., 2014).

Fukui Functions
Evaluation of the Fukui functions has been employed to explore the local reactivity of the molecules. Yang and Mortier (1986) defined the Fukui function as the first derivative of the electronic density ρ(r) of a system with respect to the number of electrons (N) at a fixed external potential ν(r), as given Equation (13). Roy et al. (1999) defined the electrophilic and nucleophilic Fukui functions for a site k in a molecule by using left and right derivatives with respect to the number of electrons as expressed via Equations (14-16).
For radical attack (16) where ρ k (N), ρ k (N-1), and ρ k (N+1) are the gross electronic populations of the site k in neutral, cationic and anionic system, respectively.
As it is known, the concept of generalized philicity have been introduced by Chattaraj et al., they defined a local quantity called philicity associated with a site k in a molecule with the assistance of corresponding condensed-to-atom variants of Fukui function, f α k as in Equation (17) (Parthasarathi et al., 2004).
where α = +, − and 0 corresponds to local philic quantities describing nucleophilic, electrophilic and radical attacks, respectively. In the light of Equation (17), the highest ω α k + corresponds to the most electrophilic site in a molecule. In addition, Lee et al. (1988) proposed different local softness, which can be used to describe the reactivity of atoms in molecules, which can be defined as in Equation (18).
Recently, Morell et al. (2006) put forward a dual descriptor, f (k), which is defined as the difference between the nucleophilic Likewise, the associated dual local softness have also been defined as expressed in Equation (20).
It is also defined as the condensed version of f k multiplied by the global softness σ . The multiphilic descriptor, ω k , is defined as the difference between the nucleophilic and electrophilic condensed philicity functions. This parameter may be used as an index of selectivity toward nucleophilic attack, which can as well characterize an electrophilic attack and is given by Equation (21) (Padmanabhan et al., 2006).
If ω k (or f (k)) > 0, the site k is favored for a nucleophilic attack, whereas if ω k (or f (k)) < 0, the site k may be favored for an electrophilic attack.

Molecular Dynamic Simulation
Adsorption characteristics of studied thiophene derivatives on iron metal surface are investigated by molecular dynamic simulation employing the Forcite module in Materials Studio 8.0 software. As model iron surface, Fe(110) surface was considered because it possesses a density packed surface and is the most stable among three common iron substrates (Guo et al., 2014). The Fe(110) system was simulated through a repeated supercell containing a 5-layer slab of Fe with 80 atoms per layer, in a 8 × 10 two-dimensional periodicity. A vacuum region 40 Å thick was included between repeated surface slabs. To simulate the metal-inhibitor systems, COMPASS force field (Sun et al., 1998) was used. The simulations of three thiophene derivatives labeled as A, B, and C on iron surface were carried out to determine the optimal adsorption sites for these molecules. All simulations made in the study were performed in an NVT canonical ensemble at 298 K with a time step of 1.0 fs and a total simulation time of 1,000 ps. The operational temperature was monitored via the Andersen thermostat. In the calculations, vacuum media was preferred and five layers of iron atoms were used.
To calculate the adsorption energies (E ads ) on Fe(110) surface of modeled of A, B, and C molecules, we used the Equation (22). Herein, it is important to note that binding energy (E binding ) is defined by the negative value of adsorption energy as given in Equation (23).
where E complex is the total energy of an inhibitor molecule and the metal surface system. E Fe is described as the energy of iron surface

Global Reactivity
As mentioned in the computation section, the investigated inhibitors were optimized by performing two different methods and using three different basis sets. The optimized electronic structures correspond to energy minima with no imaginary frequencies. According to frontier molecular orbital theory and in consistent with Fukui's theory, a high E HOMO value means the ability of a molecule to donate electrons to an assigned acceptor (metal surface in our case) with empty molecular orbital that facilitated the adsorption process and therefore indicated good inhibition performance (Obot et al., 2015). In contrast, E LUMO is related to electron affinity, which corresponds to a tendency for electron acceptance. Accordingly, the gap between energy levels of the molecules ( E = E LUMO -E HOMO ) is a significant descriptor that need be calculated. It demonstrates inveterate electron donating ability and measures the interaction between inhibitor molecules and substrate surface. The optimized molecular structures, HOMO, LUMO, as well as the molecular electrostatic potential of the investigated molecules are showed in Figure 2.
Many investigations (Khalil, 2003;Gece, 2013;Frau and Glossman-Mitnik, 2017) reported that the inhibition efficiency of inhibitors has been found to correlate with some other quantum chemical parameters like η, σ , χ, and dipole moment (DM). Besides, ω, ε, N, E b-d , ψ, which have been calculated in terms of η and χ, are also highly useful descriptors in corrosion inhibition studies of organic molecules. Our computed quantum chemical parameters for three inhibitor molecules in gas phase and water phase are given in Tables 1, 2.
As given in Table 1, it is notable that inhibitor A has the highest E HOMO among all the studied inhibitors. This reflects that the electron-donating ability of Inhibitor A is strong. As is known, corrosion inhibitors with low E values provide better inhibition performances. This is because that the excitation energy to remove an electron from the last occupied orbital will be low. It was also reported that a molecule with a low energy gap could be more polarizable, which is usually associated with a high chemical reactivity and low kinetic stability, termed as soft molecule (Madkour and Elroby, 2015). Jafari et al. (2013) also pointed that adsorption of inhibitor molecule onto a metallic surface occurs at the site of the molecule which has the greatest softness and lowest hardness. Our results in Tables 1, 2 show that all the elected levels inhibitor A has the lowest E energy in both gas and aqueous phases, and hence the molecule could have a better inhibitive performance on back-donation), N (the fraction of electrons transferred), ψ (the initial molecule-metal interaction energy), and dipole moment (DM) values for studied thiophene derivatives in gas phase and aqueous phase.

Gas phase
Aqueous phase the iron surface as corrosion inhibitor. Based on the above discussion, we can write the corrosion inhibition efficiency (in both gas and aqueous phases) order as: A > B > C. These results are in good agreement with the available experimental results (Fouda et al., 2014). However, from the results obtained for E LUMO in gas and aqueous phases the trend is irregular, which does not correlate well in the experimentally determination inhibition efficiency. So we claimed that LUMO energies of molecules may fail in terms of the explanation of their inhibition efficiencies. Generally, chemical hardness (η) is the resistance against electron cloud polarization or deformation of chemical species. Thus, the η value of a molecule and its inhibition efficiency are inversely proportional to each other because a hard molecule is renitent to give electrons . Lukovits et al. (2001) reported that η, σ , and E are quantum chemical descriptors closely correlated with each other. As it is mentioned in computation part and according to Koopmans's theorem, both softness and hardness are obtained on the basis of HOMO and LUMO orbital energies. Hard molecules with high E cannot act as good corrosion inhibitor. Nevertheless, soft molecules which have low E could be excellent corrosion inhibitors since they can easily donate electrons to metals. Based on the results reported in Tables 1, 2, it is clear that the sequence of inhibitive efficiency for studied molecules based on their hardness, and softness values can be written as: A > B > C.
In the light of the simple charge transfer theory for donation and back-donation of charges proposed by Gomez et al. (2006), the electronic back-donation process can probably affect the interaction between the inhibitors and substrate surface. As given in Equation (11), when the electron transfer and back-donation processes occur simultaneously, the energy change is directly proportional to the hardness of the inhibitor molecule. The E b−d indicates that since η > 0, then E b−d < 0, and the charge transfer from a molecule, followed by a back-donation to the molecule, is energetically favored (Bedair, 2016). Based on this principle, it is available to compare the stabilization among the inhibitor molecules, because there will be an interaction with the same metal, it is obvious that E b−d will decrease with the hardness increases. According to our results given in Table 3, as expected and in agreement with the experimental results (Fouda et al., 2014), the calculated E b−d exhibit the tendency: Electrophilicity index (ω) represents the propensity of a molecule to receive electrons. Conversely, nucleophilicity (ε) indicates the tendency to donate or share electrons with others, it is defined as the inverse of electrophilicity (1/ω). It is generally assumed that a molecule that has a large electrophilicity value is ineffective against corrosion while a molecule that has a large nucleophilicity value is expected to an excellent candidate as corrosion inhibitor. Based on Tables 1, 2, it is evident that the inhibiors have low electrophilicity index values and are good nucleophiles. But there exists discrepancy in trend between the HF and B3LYP functionals for electrophilicity, which is due to the quadratic dependence on the electronegativity. In view of this, their molecular reactivity cannot be predicted accurately, it is necessary to take into account another additional criterion to determine their inhibitive capacity.
In this work, the number of electrons transferred ( N) of the between metal substrate and inhibitor molecules was calculated using Equation (10). The results are also gathered in Table 3. Based on the Sanderson's electronegativity equalization principle (Sanderson, 1983), the charge transfer process between metal and inhibitor will continue until their electronegativity values are equal with each other. In fact, N can be regarded as a derived descriptor from the electronegativity/hardness equalization principle. It was pointed out that the positive value of electrons transferred ( N) indicates that the molecules act as electron donors. Based on the tabulated results (Table 3), we can see that molecule A has the largest number of the fraction of transferred electrons ( N) in both gas and aqueous phases, turn is, A > B > C, regardless of phase and the elected levels. These results agree well with the experimental results. Another important parameter is the initial molecule-metal interaction energy ( ψ), which has been introduced in Kokalj's work (Kovacevic and Kokalj, 2011). In our work, we have calculated ( ψ) for all the researched inhibitors and the results are given in Table 3. The results indicate that the trend of ψ is also A > B > C.
Recently, some authors used the dipole moment as an indicator of corrosion inhibition efficiencies of molecules (Gece, 2008;Zarrouk et al., 2014). Several authors showed that corrosion inhibition efficiency increases with the increase of dipole moment (Stoyanova et al., 2002;Sahin et al., 2008). Considering the idea that increasing value of dipole moment facilitates the electron transport process. Yet others proposed the opposite correlation, that is, a low value of dipole moment favors the accumulation of inhibitor molecules on the metal surface and ultimately increasing the inhibition performances (Khalil, 2003;Lebrini et al., 2005). As shown in Table 3, the calculated the dipole moments for studied compounds are irregular. There is no any remarkable relationship between dipole moment and inhibition efficiency. Thus it cannot be used a priori to judge the inhibition effectiveness.
Overall, based on the global descriptors considered for each molecule shown in Tables 1-3, The order of chemical activities from high to low is: -C(=O)-NH 2 > FIGURE 4 | (A) Score and (B) dendrogram plots obtained for three studied inhibitors.
-C(=O)-OC 2 H 5 > -C≡N substituents. Correspondingly, the inhibitive effectiveness order for the thiophene molecules is: InhA > InhB > InhC. This indicates that our calculated theoretical results are in agreement with experimental orders.

Local Reactivity
In order to have an understanding on the local reactivity of the thiophene derivatives, the Fukui indices for every atom in the inhibitors have been calculated at the B3LYP/6-31++G level. It is well known that an analysis of the Fukui indices and the local descriptors provides a more comprehensive information of the reactivity of the molecules under probe. To complete the picture, the local softness, local electrophilicity, and the dual descriptors have been also calculated for each atoms in the studied molecules. Figure 3 represents graphically the local dual descriptors, f k , σ k , and ω k for the three compounds. It should be noted that the numbering of atoms, which is given in Figure 1 are employed in this analysis. Generally, the condensed Fukui functions can make us to distinguish each part of the inhibitor molecule in the light of its distinct chemical behavior with different substituent functional groups. Therefore, the site for nucleophilic attack will be the place where the value of f − is a maximum. Conversely, the site for electrophilic attack is controlled by the value of f + .
The results reported in Tables S1-S3 reveal that, for nucleophilic attack, the highest f − k values of InhA are N19, N17, O20, and C16 atoms. For InhB, the most nucleophilic sites are N17, C16, and C22 atoms. For InhC, the most reactive sites are C16 and N17. This indicates the propensity to donate electrons to vacant molecular orbital on the iron surface to form coordinate bond. This agrees with the results of the computed HOMO density. For electrophilic attack, the highest f − k values of the three studied inhibitors are S1, C18, C15, and C13, indicating that the sites most capable for an electrophilic attack that is through which the molecule accepts electrons to form feedback bonds with Fe(110) surface. This also conforms to the computed LUMO orbital density. These results are also supported by the values of the local dual indices ( f, σ , and ω), which indicate that these inhibitors have many active sites and most of these centers have values of the three descriptors of lower than 0, except some atoms, which found to be have values > 0 (see Figure 3), indicating an electrophilic centers. A close inspection would reveal that all the molecules had the back-donation process at their carbon atoms in agreement with the frontier orbital results obtained. According to these results, one can conclude that Inh A molecule will have many active centers to interact with iron substrate. These are most likely those areas that containing N and O atoms, which are the most possible sites for bonding to iron surface through donating electrons to the Fe 3d orbitals (Khaled, 2010). Also, it can be suggested that the binding between the surface of the metal with the InhA is stronger than that in the case of B and C, respectively. Finally, the above local descriptors reveal that the theoretical order for the variation of inhibition efficiencies of the investigative inhibitors agrees with the available experimental data and it is as follows: A > B > C.

PCA and AHCA Analysis
In this work, all calculated variables have been auto scaled to compare them at the same level. Thereafter, principal component analysis (PCA) was adopted to reduce the number of variables and select the most relevant ones, which are responsible for the reactivity of researched thiophene derivatives. After performing many tests, a good separation was obtained between more active and less active thiophene compounds using 11 variables: E HOMO , I, E, χ, µ, ï, σ , ω, N, ψ, and E b−d . As indicated from PCA results, the first two principal components (F1 and F2) describe all of the overall variance as follows: F1 = 84.58% and F2 = 15.42%. The score plot of the variances is a reliable representation of the spatial distribution of the points for the data set studied after explaining almost all of the variances by the first two principal components.
In Figure 4A, the most informative score plot for the inhibitors is presented (F1 vs. F2). It is evident from the figure that PCA is responsible for the separation between more active InhA and InhB and less active Inh C where F1 > 0 for the more active compounds, and F1 < 0 for the less active one. These results are in a well agreement with the experimental results, the calculated global and local descriptors. Figure 4B shows AHCA analysis for the inhibitors under probe. The horizontal lines represent the inhibitors and the vertical lines the similarity values between pairs of inhibitors, an inhibitor and a group of inhibitors and among groups of inhibitors. It is noticed that AHCA results are very similar to those obtained with the PCA analysis, i.e., the studied inhibitors were grouped into two categories: More active inhibitors A and B and less active one C.

Molecular Dynamic Simulations
Molecular dynamic simulation approach is very important in terms of the explanation of the nature of the interactions between corrosion inhibitors and metallic surface. The optimized equilibrium adsorption configurations for A, B, and C molecules on Fe(110) surface are given in Figure 5. Adsorption energy is known as the energy released when inhibitor molecule was adsorbed on metal surface. As given in previous section, the binding energy is the negative value of the adsorption energy. Higher negative value of adsorption energy and higher positive values of binding energy represent the more stable and more strong interaction between metal surface and inhibitor molecule. In Table 4, calculated adsorption and binding energies as well as experimentally determined corrosion inhibition efficiencies for studied thiophene derivatives are presented. It is apparent that the binding energies of three derivatives on the Fe(110) substrate decrease in the order A > B > C, which is in consonance with the experimental inhibition efficiency orders (Fouda et al., 2014).

Comparison Between Experimental and Theoretical Results
In this subsection, a comparison between experimental and theoretical results is presented. Our study shows that there is an excellent correlation between our theoretical results (global quantum descriptors and MDs results) with the experimental inhibition efficiency (IE%). Figures 6A,B shows graphical representation of the relationship between the linear correlation R obtained for the relationship between theoretical reactivity parameters calculated (in both gas and aqueous phases) using B3LYP and HF methods with 6-31++G basis sets for the studied inhibitors and their experimental IE%. The results of linear relationship coefficients are summarized in Table S4. However, the results obtained via B3LYP method are well correlated with the experimental results than those obtained by HF method. As can be seen in Figure 6, there are a very excellent linear correlation between the experimental inhibition efficiency and the theoretical descriptors in gas phase for both methods in gas phase. It is apparent from the column graphs plotted that B3LYP method provided more accurate results compared to HF method. In addition, it can be said that calculation levels including 6-31++G basis set is more successful compared to other calculation levels in terms of the obtaining good agreement with experimental results. Finally, it can also be seen from Figure 7 that there exists closely correlation between experimental anticorrosion efficiencies and calculated binding energies with a high correlation coefficient of 1.00.

CONCLUSIONS
In this work, Hartree Fock as well as DFT calculations, molecular dynamic simulation, PCA, and AHCA were used to analyze the anticorrosive performances of some thiophene derivatives against the iron metal. Global and local reactivity descriptors were calculated in both gas and aqueous phase for three studied inhibitors. Within the framework of the theoretical results obtained in this study, the following conclusions can be written.
(1) DFT, molecular dynamic simulation, PCA, and AHCA results showed that corrosion inhibition efficiency ranking of studied molecules is given as: InhA > InhB > InhC.
(2) It is apparent from binding energies and adsorption energies calculated for studied thiophene derivatives, these molecules are very effective against the corrosion of iron. (3) According the PCA and AHCA results, least active inhibitor among the studied molecules is inhibitor C. (4) Both theoretical data and experimental results are compatible with inductive effect of functional groups appearing in the molecular structures of studied thiophene derivatives. (5) Theoretical results obtained in this work have far-reaching significance to the rational design of novel thiophene derivatives as corrosion inhibitor. Frontiers in Chemistry | www.frontiersin.org