Investigation of Corrosion Inhibitors Adsorption on Metals Using Density Functional Theory and Molecular Dynamics Simulation

The use of computational chemistry as a tool in the design and development of organic corrosion inhibitors has been greatly enhanced by the development of density functional theory (DFT) and Molecular dynamic simulation (MD). Experi-mentally corrosion inhibitor development requires lots of money and time. Thus, in the era of hardware and software development, corrosion scientist can select a potential inhibitor on the basis of theoretical analysis of molecular properties of inhibitor molecules, which have reduced the cost. DFT and MD are capable to accurately predict the inhibition characteristics of inhibitor molecules using molec-ular/electronic properties and reactivity indices. The purpose of this book chapter is to summarize some important features related to DFT and MD, giving a brief background to the selected DFT/MD-based chemical reactivity concepts, calculations and their applications to organic corrosion inhibitor design. The impact of this book chapter is to illustrate the enormous power of DFT and MD.


Introduction
Corrosion of active metals such as iron (Fe), copper (Cu) and Aluminum (Al) can be reduced by modifying their surface using organic/polymeric corrosion inhibitors. The inhibitor molecules adsorbed onto the metal surface and form a protective barrier against corrosion [1]. The corrosion protection ability of inhibitors depends upon the extent of interaction between inhibitor molecules and metal surface [2]. The basics mechanism through which adsorbed inhibitor molecules alter the corrosion process consists of either changing the cathodic and anodic corrosion reactions or physically blocking the active sites present on metal surface or by [3].
Perusal of the literature clearly reveals that huge number of publication exits where experimental works have been carried out to understand the corrosion inhibition process, however use of theoretical studies including density functional theory (DFT) and molecular dynamic simulation (MD) appeared in recent years. In traditional experimental work scientists has to test huge number of organic compounds for selecting them as a potential corrosion inhibitor. This kind of search takes more time and money. However, in the era of development in the field of hardware and software technology make us capable to select organic compounds as potential corrosion inhibitor without expensing huge amount of money on chemicals. The corrosion inhibition efficiency of organic compounds based on molecular and electronic properties can be predicted accurately using MD and DFT calculations. In the recent years, scientist have devoted their positive efforts in the field of corrosion using theoretical calculation [4][5][6]. To our knowledge, the review exists in the literature only describing separately DFT and MD and not both in relation to corrosion inhibition study [7,8]. Here we are trying to describe both DFT and MD on the same platform simultaneously.
In this chapter we will try to explore some recent studies and concept which can provide some new concept in the field of computational chemistry related to corrosion inhibition phenomenon. Here, we summarize the adsorption of organic molecules onto the metal surface especially iron, aluminum and copper in order to study their corrosion inhibition properties using DFT and MD simulations.

Basic concepts related to molecular modeling
In literature, large number of publication/papers exits where corrosion inhibition property of inhibitors has been studied using theoretical study. In general organic molecules have the ability to donate electrons to empty d-orbital of metal surface and form covalent bond which is known as coordinate bond. However, they also undergo back donation, where they accept the free electrons presenting in filled metal orbitals and this is called as retrodonating bonds. Thus, these kinds of bonds make the organic compounds a potential candidate for corrosion inhibition. The theoretical parameters like highest occupied molecular orbital (HOMO), lowest unoccupied molecular orbital (LUMO), dipole moment (μ), energy gap (ΔE), Fukui indices, charge density, polarizability, softness, etc. [9-11] were discussed in DFT calculation.
DFT study separately represents the organic compound and metal surface [12][13][14][15][16]. In DFT, the interaction between the inhibitor molecule and metal surface was not directly modeled. However, in MD the most reactive fragment of the molecule that has greater affinity to interact with the metal surface is study using ab initio method. This method provide a real picture which is actually going in the corrosion process like interaction between inhibitor molecule and metal surface, orientation of molecule towards metal and organization of molecules.

The basics of DFT: the Hohenberg-Kohn theorem
Quantum electrochemistry is the field which includes quantum mechanics, electrochemistry and electrodynamics. In general, quantum electrochemistry is an application of density functional theory (DFT) for studying the electrochemical processes like transfer of electron towards electrode surfaces [17].

Functional
The functional is a function of a function, i.e., electron. Therefore, on the basis of fundamental quantum mechanics and from the experimental results different functional have been developed and are listed in Table 1 [29].
There are several generations of density functional. First generation is known as X α method and it is the simplest one. This functional was developed by J.C. Slater, who was working on the approximation of Hartree-Fock but unwillingly discovered the simplest for DFT. In X α functional electron exchange but not the correlation was included. The results obtained from X α method is as accurate as HF but in some cases it is better than HF.
The functional of second generation uses both the density and its gradients. In 1986 first gradient-corrected energy functional was proposed in which Becke, Perdew and Wang proposed exchange and Perdew correlation [32,33]. In the present time the most popular functionalis of Becke and Perdew for exchange and correlation respectively [32], also Lee et al. [34] for correlation and Perdew and Wang for exchange and correlation both [35]. These functional are collectively known as "generalized gradient approximations (GGAs)." The third generation is hybrid functional and it is more advanced than GGA. In this, the functional consists of both Hartree-Fock type exchange and DFT exchange calculated from the orbitals. In 1998, a hybrid functional called as B3LYP functional was introduced. In 1994, Gaussian introduced this functional in in computational package for first time and it is written by below equation: where E X(HF) represent the Hartree-Fock exchange energy, E X(LSD) stands for Dirac exchange energy and E X(GGA) and E C are the gradient corrections to exchange energy, i.e., Be88form and Lee-Yang-Parr (LYP) functional respectively [34,36].

Basis sets
In DFT calculation various basis sets are used [37]. The simplest one is STO-3G. In this basis set, 1s is given by three Gaussians and 2s, 2px, 2py, and 2pz each by another three. Improvement in the basis set can obtain by incorporating two 1s functions for hydrogen and for 2nd row atoms like carbon two 2s and two 2p functions has incorporated. These functional are called as split-valence basis sets.
Here valence and core orbitals are represented by two sets of functions and single set of functions respectively. An example for this: Carbon 3-21G: Three Gaussians for 1s combination, same two Gaussians for 2s and 2p combination; plus 2s' and 2p' the same one Gaussian. Carbon 6-31G: Six Gaussians for 1s combination, same three Gaussians for 2s and 2p combination; plus 2s' and 2p' the same one Gaussian.
In general it can be represented as i-jk, here i denote the number of Gaussians that are representing each core basis function, and j and k represent the numbers of Gaussians for split-valence basis functions.
Further advancement generates triple-split-valence basis sets, for example 6-311G. In this basis set core functional consists of six Gaussians and three sets of valence functions containing three, one, and one Gaussians, respectively.
Multiple zeta basis sets such as double zeta (DZ) and triple zeta (TZ) basis sets has also introduced. The difference between multiple and split-valence basis sets is that in multiple basis sets all of the orbitals are separated into either two or three sets of functions, and so forth. Also, multiple basis sets used different α-coefficients for s and p orbitals.
So, improvement in the electronic structure calculation can be achieved by adding functions corresponding to orbitals with a higher angular momentum. This concept

Basis
Options Atoms can explained as follows: for hydrogen atom, p functions are added; for carbon, nitrogen, oxygen, etc., d functions are added; and for transition metals, f functions are introduced. These are represented by adding an asterisk to the basis set and also specifying the p, d, f functions. This can be understood as per below example: 6-31G * or 6-31G(d): Adds d functions to 2nd row elements (C, N, O, etc.). 6.31G * * or 6-31G(d,p): Adds d function to 2nd row elements (C, N, O, etc.) and p functions to H.
The lone pair of electrons on the heteroatoms can be taken into confederation by adding diffuse functions and these are represented in the basis set using + and ++ signs. For hydrogen atom polarization and diffuse functions are not very import but for other atoms like carbon, nitrogen, oxygen it is important, for example, 6-31+G: represents that diffuse functions has been added to 2nd row elements (C, N, O, etc.). 6-31++G: represents that diffuse functions has added to 2nd row elements and H.
In the case of transition elements which contain large number of electrons, requires more time for calculation. So, the LanL2DZ basis set is used for transition metal calculations. This basis set is derived by combining the valence electrons using double zeta functions and Los Alamos ECP. The software used in corrosion inhibition study with a wide range of basis sets are Spartan, Material Studio, Gaussian 03, Gaussian 09 and most recently Gaussian 14 has been introduced.

DFT parameters and their application in corrosion inhibition study 4.1 Frontier molecular orbitals
The optical and electronic property of organic compounds is explained by analyzing their frontier molecular orbitals (FMO) like highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) [38]. The interaction of organic compounds as corrosion inhibitor with the metal surface depends upon the electron donation ability of the organic compounds, i.e., ionization potential (I.P.). However, the electron accepting ability of organic compounds corresponds to the electron affinity (E.A.) [39]. So, HOMO is the electron donating orbitals and LUMO is electron accepting orbitals [40]. In general for a good corrosion inhibitor HOMO value should high, LUMO should low and energy gap between the LUMO and HOMO (ΔE) should low [40].
The HOMO and LUMO are represented in Figures 2 and 3. From the figure, it could be said that HOMO regions are distributed on the triazine ring, phenyl ring containing electron donating group, nitrogen and oxygen atoms. These regions correspond to the electron donating areas and inhibitor molecules interact with the metal surface through these regions and undergo adsorption over it. However, the LUMO are distributed over the triazine and phenyl ring respectively and these are corresponding to the electron accepting regions. From these LUMO regions inhibitor accept electron from metal. Order of inhibition efficiency is HT-1 > HT-2 > HT-3.

Electronegativity and the electronic chemical potential
One of the important global molecular property calculated by density functional theory (DFT) which provide an idea for chemical reactivity and selectivity of organic compound as corrosion inhibitor are electronegativity (χ) and chemical potential (μ) [46]. Mathematical form of this is as follows [47,48]:   where I.P. and E.A. are ionization potential and electron affinity.

Global hardness and softness
According to Parr et al. global hardness (η) measures the resistance of an atom to transfer the charge [49] and it can be calculated using below equation: where E L and E H are the HOMO and LUMO energy. The reciprocal of global hardness is global softness (σ) and given as [50]:

Electric dipole polarizability
The electric dipole polarizability (α) is a measurement of the linear response of the electron density in the presence of an infinitesimal electric field F and it describes a second-order variation in energy: The polarizability (α) is calculated as the mean value and express through above equation: It was discovered that polarizabilities are inversely proportional to the third power of the hardness values [51][52][53][54].

Electrophilicity index
Parr et al. has introduced the concept of global electrophilicity index (ω) and is given as follows [55,56]: The global electrophilicity index (ω) is also related to the ionization potential (I.P.) and electron affinity (E.A.) and is given as follows [56]: The electrophilicity index accounts the ability of the molecule to accept electrons. So, higher the ω value, greater would be the molecule to accept the electrons.

The fraction of electrons transferred
When the inhibitor and metal interaction occurs, flow of electron takes place from the lower electronegativity molecule to the higher electronegativity metal. This transfer of electron continues until the chemical potential becomes equal. The fraction of electron transfer is given as follows [57]: where metal and inhibitor molecule are represented by m and i indices respectively. The absolute electronegativity of metal and inhibitor molecule are represented by χ m and χ i respectively. The absolute hardness of metals and inhibitor molecule is given by ηm and ηi, respectively. In case of iron the theoretical values of χ m and ηm are 7 and 0 eV mol À1 respectively [57]. However, the value of 7.0 eV corresponds to the free electron gas Fermi energy of iron in the free electron gas model. But the use of this value for χ m is conceptually wrong because here electronelectron interaction has been neglected [58]. Thus, work function (φ m ) of the metal has been incorporated for an appropriate measurement of electronegativity. So, the new formula for the estimation of N is as follows [59]: The work function (φ) values for iron surface plan Fe (1 0 0), Fe (1 1 0) and Fe (1 1 1) surfaces are 3.91, 4.82 and 3.88 eV, respectively [59,60]. The surface energy of iron planes for bcc structure are in order of Fe (110) < Fe (100) < Fe (111) and out of these the most stable on plan is Fe (110). Thus, most commonly Fe (110) plan is used in corrosion inhibition research [61].

Fukui function and local softness
In the field of corrosion research inhibition of corrosion is governed through donation and acceptance of electron, which basically involves the nucleophilic and electrophilic reactions. This can be achieved through the evaluation of the Fukui indices [61][62][63]. Thus, by invoking the HSAB principle in a local sense, one may establish the behavior of the different sites with respect to hard or soft reagents.
The Fukui function generated using the finite difference approximation is as follows: where f + is nucleophilic and f À is electrophilic Fukui functions respectively, ρ (N+1) , ρ (N) and ρ (N-1) are the electronic densities of anionic, neutral and cationic forms of the atoms with N+1, N and N-1 electrons.
Castro et al. [64] studied the series of isomers of the Schiff base (E)-2-(2-hydroxybenzylideneamino) phenylarsonic acid as corrosion inhibitor. The author has pictorially presented the negative Fukui, i.e., f À function and is given as follows.
In Figure 4, the black region represents the area where inhibitor contains the major portion of electrons for donating to the empty d orbitals of Fe metal. The Fukui function represents the change of the electronic density in a given point with respect to the change in the number of electrons N [65]. For the electrophilic attack to occur, the contributing atoms are likely to be the carbon atoms on the OHsubstituted phenyl ring in all isomers, the oxygen atom in the isomer SBAs-OH-I, the nitrogen atom of the imine group of isomer SBAs-OH-III, and a region located in AsO(OH)2-substituted phenyl ring.
Local softness is calculated is using the equation [66]: The integration of equation gives global softness value: Local softness and Fukui function are related as per below equation: The local softness contains information similar to those obtained from Fukui function plus additional in-formation about the molecular softness, which is related to the global reactivity with respect to a reaction partner.

Charge transfer: donation and back donation
The energy change in the molecule is attributed because of two processes either electron donation to the metal or acceptance from the metal, i.e., back donation. Gomez et al. [67] proposed a simple charge transfer model for donation and back donation. In this model when the molecule receives certain amount of charge (ΔN À ) the change in energy is given as follows: However, when the molecule back donate certain amount of charge (ΔN À ) the change in energy is given as follows: If the total change in energy (ΔE T ) is the summation of ΔE + and ΔE À , and assuming that the donation and back donation, i.e., (ΔN À ) and (ΔN + ) are equal, than: In a situation when the total energy change becomes a minimum with respect to ΔN + , then and change in total energy becomes as follows: According to Guo et al. [68] the change in total energy are negative, i.e., Δ E T < 0. This implies that the back donation of charge from metal surface to inhibitor molecules is favorable.

Inhibitor/surface interaction study using molecular dynamic simulation (MD)
Molecular dynamics is a computer-based modeling technique by which the evolution as a function of time or trajectory of a molecule is described by the principles of classical Newtonian mechanics [69]. In molecular dynamics, intramolecular movements of inhibitor molecules are simulated in order to visualize the real picture of corrosion inhibition process. Molecular dynamics simulation calculation has two serve two purposes: a. Intramolecular movements' simulation and calculation of thermodynamic properties such as entropy, free energy, etc.
b.Optimization of the molecular structures efficiently by avoiding multiple minima by annealing the simulation process. The range of simulation time is approximately between 10 À14 and 10 À10 s per simulation.
Molecular dynamics simulation involves the implementation of initial steps: In the first very step build molecular structure are introduced or immersed in a box having solvent molecules, after that simulation process starts. The total energy of the molecular system should be minimizing before the start of the simulation, in order to avoid the generation of aberrant trajectories (when the initial forces are too great). For achieving this, the steepest slope/steepest descent method and that of the conjugated gradients have been frequently used.

Adsorption behavior of inhibitors on mild steel surface
In literature, several DFT calculations are observed where inhibition efficiency and molecular structure/electronic properties of organic corrosion inhibitors have been correlated. Xu et al. [70] has studied corrosion inhibition property of Schiff's base derivatives via DFT calculation using GGA/BLYP method. The result indicates that 2-PCT molecules adsorb strongly onto iron surface and provide higher inhibition efficiency.
Zhang et al. [71] investigated the inhibition of methionine (Met) and proline (Pro) in PCMs solution using DFT and B3LYP/6-31G (d) method. The inhibition efficiency of Met and Pro follow the order of Met > Pro.
The corrosion inhibition effects of thiourea (TU), methylthiourea (MTU) and phenylthiourea (PTU) on mild steel in 0.1MH2SO4 was investigated by Ozcan et al. [72]. The electronic properties were calculated using ab initio RHF/6-31G(d) method. The energy of E HOMO , the energy of E LUMO and energy gap values supports that PTU is best inhibitor.
The corrosion inhibition performance of three naphthyridines derivatives were explored for mild steel in 1 M HCl by Singh et al. [73] using DFT methods. The inhibition efficiencies of N-3 and N-2 are higher due to the presence of electrondonating -OCH 3 and -CH 3 substituents respectively as compared to N-1, which has no substituents.

Adsorption behavior of inhibitors on aluminum surface
The inhibitive action of poly ethylene glycol (PEG) against the corrosion of aluminum surface in 1 M HCl was investigated by Awad et al. [75] using B3LYP/6-31 + G(d,p) basis set. DFT was performed on repeating units of 1, 2, 3, 4 and 5. The inhibition performance increases with increase in number of repeating unit.
The author Li et al. [76] studied the corrosion inhibition effect of three oxime compounds on the corrosion of aluminum in HCl solution using BLYP in conjunction with double numerical plus d-functions (DND) basis set. In these inhibitors, HOMO and LUMO electron densities are localized principally on C=N-OH, which suggest that C=N-OH functional group could both accept and donate the electron.
Khaled et al. [77], investigated the corrosion effect of and adsorption characteristics of three imidazole derivatives on aluminum in 1.0 M HCl using local density functional (LDF) method with a double numeric polarization (DNP) basis set, and a Becke-Perdew (BP) functional. A good correlation between the rate of corrosion and EHOMO, as well as with energy gap (ΔE = E LUMO À E HOMO ).
Kaya et al. [78] applied HF and DFT/B3LYP methods with SDD, 6-31++G (d, p) and 6-31 G basis sets both in gaseous and aqueous phase for studying the adsorption and corrosion inhibition properties of some benzotriazole and phospono derivatives on aluminum. The inhibition efficiency ranking of these molecules as: PBA > PBTA>PAA > TBTA.

Adsorption behavior of inhibitors on copper surface
Corrosion inhibition of indazole (IA) and 5-aminoindazole (AIA) for copper in NaCl solution was studied by Qiang et al. [79]. The molecular structures of IA and AIA were geometrically optimized by density functional theory (DFT) using B3LYP functional with 6-311++G(d,p) basis set in aqueous phase. The higher value of HOMO and lower value of ΔE in case of AIA concludes that the interaction of AIA with copper will be stronger than that between IA and copper. Khalid and Madkour et al. [13,80,81] probe the interactions between benzotriazole, methionine, phenol derivatives and copper via quantum calculations respectively. The obtained results suggests that quantum parameters are very useful in characterizing organic compounds as an adsorbate.

Corrosion inhibitors studied using MD
Feng et al. [74] performed the MD simulation of 1-[N,N 0 -bis (hydroxylethylether)-aminoethyl]-2-stearicimidazoline (HASI). The MD calculation suggests that adsorbed imidazoline molecule is parallel to the iron surface in order to maximize its contact with the surface (Figure 5). The interaction energy for Fe atom is À284 kJ/mol, for Fe 3 O 4 is À226 kJ/mol and for Fe 2 O 3 is À157 kJ/mol. So, inhibitor adsorbed on Fe surface more strongly than Fe 3 O 4 and Fe 2 O 3 .
MD calculation of adsorption of indazole (IA) and 5-aminoindazole (AIA) on copper surface was studied Qiang et al. [79]. The optimized equilibrium configuration of inhibitors molecules are shown in Figure 6 and inhibitors adsorption occurs through parallel orientation. The adsorption energy between Cu (1 1 1) surface and inhibitor molecule are À250.86 kJ/mol for IA and À307.86 kJ/mol for AIA.
Molecular level MD calculation of Benzotriazole and Phospono derivatives was studied by Kaya et al. [78] on aluminum surface. The best adsorption configuration of inhibitor molecules on Al (111) surface are represented in Figure 7. The parallel configuration of inhibitor molecules suggests the stronger adsorption.

Conclusion and future work
Computational program provide an atomistic view for understanding the corrosion inhibition process and significantly contribute for the selection of adequate inhibitor molecules. DFT calculations has made possible to characterize a suitable inhibitor molecules using their intrinsic properties such as energies of HOMO and LUMO orbitals, energy gap (ΔE), electronegativity, hardness, softness, etc. It helps corrosion scientist to understand the inhibitor-metal interaction through the locations of HOMO and LUMO orbitals. The examination of local reactivity reveals the nucleophilic and electrophilic centers which react with the metal surfaces.
The literature survey suggests that with help of quantum chemistry a good understanding of corrosion system can be achieved when various metal surfaces and various inhibitor molecules are compared. DFT has ability to perform  calculations for complex or large organic molecules that can be applied as a potential corrosion inhibitor. In addition, DFT provide a platform to scientist for designing a unique inhibitor molecule and understand their molecular structure at atomistic level.
In recent years, computational modeling like DFT and MD provide real inhibitor-metal interaction picture by performing calculation in aqueous phase. Such an approach act as a benchmark for experimental study in the field of corrosion research.
Molecular dynamic simulation (MD) gives a representative image of inhibitor molecule orientation on the metal surface. This helps to analyze the inhibition activity of inhibitor molecules, i.e., inhibitor with parallel or flat orientation with respect to metal surface will cover more surface area and provide a better corrosion inhibition protection as compared to the inhibitor which have vertical orientation.
Future research includes modeling the electrochemical potential, generation of alloy surface layers, addition of thin oxide film covering onto the metal surface, calculation of more complex molecule, comprehension of solvent molecule role in adsorption, understanding the defect role in adsorption process, etc.