MOLECULAR DYNAMICS SIMULATIONS AND QUANTUM CHEMICAL CALCULATIONS FOR THE ADSORPTION OF SOME IMIDAZOLINE DERIVATIVES ON IRON SURFACE

The molecular dynamic (MD) simulation and quantum chemical calculations for the adsorption of [2-(2-Henicos-10enyl-4,5-dihydro-imidazol-1-yl)-ethyl]-methylamine (HDM) and 2-(2-Henicos-10-enyl-4,5-dihydro-imidazol-1-yl)-ethanol (HDE) on iron surface was studied using Materials Studio software. Molecular dynamic simulation results indicate that the imidazoline derivative molecules uses the imidazoline ring to effectively adsorb on the surface of iron, with the alkyl hydrophobic tail forming an n shape (canopy like covering) at geometry optimization and at 353 K. The n shape canopy like covering to a large extent may prevent water from coming in close contact with the Fe surface. The quantum chemical calculation based on the natural atomic charge, the frontier molecular orbital and the Fukui indices values and plots shows the active sites of the molecules to be mainly the N=C-N region in the imidazoline ring, others include the nitrogen and oxygen heteroatoms in the pendant part and the double bonded carbon atoms in the hydrophobic tail of the imidazoline derivative molecules. The quantum chemical calculations also reveal that the amine group in HDM and the hydroxyl group in HDE which is attached to the imidazoline ring do not result in a significant increase in the HOMO nor the LUMO density which can aid adsorption.HDM has a lower energy gap of 4.434 eV and 3.824 eV, a higher EHOMO of -4.273 eV and -4.152 eV and a higher global softness of 0.45 and 0.52 compared to HDE which have an energy gap of 4.476 eV and 4.084 eV, a EHOMO of -4.349 eV and -4.607 eV and a global softness of 0.45 and 0.49 at geometry optimization and at 353 K. The adsorption ability of the molecule is given as at geometry optimization HDM > HDE and at 353 K HDM > HDE. Theoretically HDM is a better inhibitor than HDE. The adsorption ability of the molecule is in line with the binding energy at the temperature studied.


INTRODUCTION
Amines and their salts are said to be the most economic and effective corrosion inhibitors for oil and gas wells (Schauhoft and Kissel, 1999). Imidazoline is a typical amine-nitrogen compound which is heterocyclic in nature. Heterocyclic organic compounds consisting of a π-system or heteroatoms such as nitrogen (N), sulphur (S), or oxygen (O) are said to be good corrosion inhibitors for metals in acidic media (Okafor and Zheng 2009, Okafor et al. 2010, Okafor et al. 2010a, Obot and Obi-Egbedi 2008. The lone electron pairs in the hetero atoms and the planarity of a molecule are important features that determine the adsorption of molecules on a metallic surface (Abd El-Rehim et al, 1999). The inhibition efficiency has been closely related to the inhibitor adsorption abilities and the molecular properties for different kinds of organic compounds (Wang et al, 2004). The power of the inhibition depends on the molecular structure of the inhibitor. Organic compounds, which can donate electrons to unoccupied orbital of a metal surface to form coordinate covalent bonds, can also accept free electrons from the metal surface by using their anti-bonding orbital to form feedback bonds, constitute excellent corrosion inhibitors (Udhayakala et al, 2012). Recently the effectiveness of an inhibitor molecule has been related to its spatial as well as electronic structure (Stoyanova and Peyerimhoff, 2008) on which Quantum chemical methods have already been proven to be very useful in determining the molecular structure as well as elucidating the electronic structure and reactivity of molecules (Kraka and Cremer, 2000) while molecular dynamic (MD) simulation has shown to give a better understanding of how inhibitor molecules behave on a metal surface. This study (MD) is very important for understanding shape changes, interactions and energetics of molecules. In this study two imidazoline derivatives consisting mainly of the 'head' (heterocyclic ring), the 'anchor' (or 'pendant') the substituent attached to the heterocyclic ring, and the 'tail' parts (the long alkyl hydrophobic tail) (Hassazadeh, 2004) has been studied. The objective is to determine their modes of adsorption on iron surface, their adsorption ability and the active sites of this molecules using Molecular dynamic simulation considering properties such as interaction/binding energy, bond length, natural atomic charge and quantum chemical parameters such as dipole moment (), energy of deformation (Đ),energy of the highest occupied molecular orbital (E HOMO ), energy of the lowest unoccupied molecular orbital (E LUMO ), energy gap (ΔE), total energy (E), electron affinity (A), ionization potential (I), global hardness (ƞ), global softness (S), and Fukui functions f(r) at the geometry optimized and at 353 K structures

The Imidazolines
The names, abbreviations, structures, molecular weights and melting points of the imidazolines used in this work are shown in Table  1  Molecular Dynamics (MD) simulations were carried out in a simulation box with periodic boundary conditions using Materials Studio 4.0 (from Accelrys Inc.) software. The box consisted of air on surface (with 4 layers of iron atoms cleaved along the (0 0 1) plane, with the uppermost layer released and the inner layer fixed), with a dimension of 20.1Å x 8.6Å x 34.4Å and a vacuum layer of height 28.1 Å. The simulation was carried out using Forcitemodule (which is a sophisticated typical molecular mechanical tool that can perform fast energy calculations as well as reliable geometry optimization of molecule as well as periodic system) with a time step of 0.1 fs and simulation time of 5ps performed at 353 K, NVT ensemble, and the force field CVFF (Consistent Valence Force Field) was used (It is mainly used for the study of structures and binding energy, though it can also predict vibrational frequencies and conformation energy reasonably well). The interaction energy was calculated as given by Xia et al, (2008) where E interaction is the interaction energy, E total is the total energy of the iron crystal with the adsorbed inhibitor molecules. E Fe and E inhibitor are the total energy of the iron crystal and the total energy of the inhibitor molecules, respectively. The binding energy is the negative value of the interaction energy.

Computational Details for Quantum Chemical Calculation
Quantum chemistry calculations were done using two modules available in the Material studio 4.0 (from Accelrys Inc.) software. These modules are Vamp (which is a highly efficient semi empirical molecular orbital package) and Dmol 3 (which is a program which uses the density functional theory (DFT) with a numerical radial function basis set to calculate the electronic properties of molecule clusters, surfaces and crystalline solid materials from the first principle). Using the Vamp module, theoretical calculations were carried out at the restricted Hartree-Fock level (RHF) using the Hamiltonian Parametric method 3 (PM3) which is based on the neglect of diatomic differential overlap (NDDO) approximation. The density functional theory (DFT) calculations was performed using BLYP (from the name Becke for the exchange part and Lee, Yang and Parr for the correlation part)functional method and the double numeric with polarization (DNP) basis set, this basis set is used because it is the best set available in Dmol 3 .The dipole moment is the measure of net molecular polarity, which is the magnitude of the charge (Q) at either end of the molecular dipole times the distance (r) Dipole moment is used to know the distribution of electrons between two bonded atoms. The energy of deformation is the energy required to change the orientation of a molecule. These were calculated using the Vamp module available in the Material studio version 4.0 software.
According to Koopmans theory the energy of the frontier molecular orbital which are the energy of the highest occupied molecular orbital (E HOMO ) and the lowest unoccupied molecular orbital (E LUMO ) is related to the ionization potential (IE) and the electron affinity (EA) according to equation (3) and (4) (Koopmans, 1933). These were determined using the Dmol 3 module present in the Material studio version 4.0 software.
Thus the value of the global hardness (ƞ) according to Pearson, operational and approximate definitions is given in equation (5)according to Parr and Pearson, (1983) The global softness (S) is the inverse of the global hardness (ƞ) as expressed in equation (6) according to Pearson, (1988) = ƞ ………………………………………. 6 The local reactivity was analyzed by determining the Fukui indices which indicate the reactive regions, in the form of the nucleophilic and electrophilic behavior of each atom in the molecule. The Fukui function f(r) is defined as the first derivative of the electronic density ρ(r) with respect to the number of electrons N at a constant external potential v(r). Thus, using a scheme of finite differences (Lee et al, 1988), we had -( ) = ( ) -( ) ……………………..... 7 for electrophilic attack and are the Fukui negative, Fukui positive, electronic densities of anionic, electronic densities of neutral species and electronic densities of cationic species, respectively.
The N corresponds to the number of electrons in the moleculeN+1 corresponds to an anion, with an electron added to the LUMO of the neutral molecule. N-1 corresponds to the cation with an electron removed from the HOMO of the neutral molecule. The Fukui indices calculations were performed using Dmol 3 module present in the Material studio version 4.0 software. All calculations in this study were done on the ground state geometry structure (geometry optimized) and the structure at 353K. The colour codes for the atoms in the molecules studied are gray for carbon, blue for nitrogen, red for oxygen and white for hydrogen.

Molecular dynamic simulation
The close contacts between the inhibitor molecules and iron surface as well as the best adsorption configurations for the compounds at geometry optimization and at the temperatures of 353 K are shown in Figure 1 (a), (b), (c) and (d).Equilibration of the system was established by the steady average values of energy as well as temperature (Xia et al, 2008).The calculated interaction, binding energy and the distance between the inhibitor molecules and the iron surface obtained from molecular dynamics simulation are shown in Table  2. From Fig 1. we can see the head group and the pendant part of the imidazoline molecule is being attached with the iron surface, while the alkyl hydrophobic tail bends to form an n shape consequently repelling the water molecules from the metal surface. Figure 1(b) clearly shows the head group of the molecule attached plainly with the iron surface.

Interaction/binding energy
Interaction energy (E interaction ) between Fe surface and inhibitor molecules can be evaluated according to the Equation (1). The values of the interaction energy for the two (2) imidazoline derivatives are shown in Table 2. From Table 2 it is observed that the two (2) imidazoline derivatives show very low interaction energies meaning that they can adsorbed on the Fe surface (Zeng et al, 2011). The more negative the interaction energy, the more it can effectively interact with the iron surface and the higher the binding energy and also the inhibition efficiency (Roger, 2006). The values of the binding energy are shown in Table 2. HDM has the highest binding energy at 353 K. Table 2 shows the distance between each of the adsorbed imidazoline molecules and the Fe surface at 353 K. From Table 2 it is seen that 3.336 Å to 3.394 Å are the distance range observed between the imidazoline derivatives and the Fe surface. Due to the high distance value (d >3.0 Å) (Roger 2006, Hellsing 2008) observed between the imidazoline molecules and the Fe surface, it may be that the adsorption is a physical adsorption (physisorption).The average metalligand bond distances, Fe-N〉lies between 1.973 Å and 2.172 Å. These are the extreme values observed for the bond length expansion between the low-spin and highspin states of a typical Fe-N bond. The distances of separation as shown in Table 2 are 3.394 Å for HDM and 3.336 for HDE. This suggests ineffective overlap of orbitals thus weak interaction between the orbitals giving rise to bonds with relatively low energy which is possible when the molecules are physically adsorbed unto the metal surface (Oliver et al 1979, Ebenso et al 2010, Bouayed et al 1998, Gece 2008, Obot and Obi-Egbedi 2010.  Fig 2 shows the optimized geometry and equilibrium structures of HDM and HDE. The changing structure of the imidazoline molecules at 353 K reveals a change in bond length, bond angle and torsional effects. This is in accordance with Xia et al. (2008). The energies of the molecules are shown in Table 3. It is observed that the energy of the inhibitor molecules at geometry optimization is lower than the energy at 353 K for each molecule. This may be because the temperature at equilibrium increases the entropy of the atoms in the molecule which leads to an increase in the energy of the system (molecule).  Bond length analysis Bond length or bond distance is the average distance between two bonded atoms in a molecule. The bond length in Armstrong (Å) at the Geometry optimized as well as at the structure at 353 K of the molecules are shown in Fig 3 (a) and (b).It is observed that changes in the structural orientation of the atoms in the molecules at 353 K brings about a change in the bond length of each of the atom in the molecule from each other. Energy of a certain value has to be supplied in order to break any bond. This means that the closer the nuclei of the bonding atoms are a greater supply of energy is needed to separate the atoms due to large force of attraction between the atoms. In other words, short bond length requires high dissociation energies to break the bond. From Fig 3, it is observed that the heteroatom's (nitrogen and oxygen) bonded to a carbon atom show a shorter bond length than the carbon to carbon bond length at geometry optimized and at equilibrium structures(HDM {N1-C2, C3-N4, N4-C5, C6-N7 and N4-C8} and for HDE {O1-C2, C3-N4, N4-C5, C6-N7 and N4-C8}) It means that more energy will be required to break this bond and hence a high chemical reactivity. It is also observed that atoms which are connected by double bonds show shorter bond length than atoms connected by single bonds (HDM {N7=C8 and C18=C19} and for HDE {N7=C8 and C18=C19}. These atoms will also need a higher energy to break its bonds and hence a high chemical reactivity. The shorter the bond length the higher the bond energy and the higher is the reactivity of the bond.

Natural atomic charge
From classical chemical theory, all chemical interactions can be done by orbital interactions or by electrostatic interaction. The presence of electric charge in the molecule is the propelling force of an electrostatic interaction. Local electric charges have been proven to be vital in several chemical reactions as well as physiochemical characteristics of compounds (Xia et al, 2008). Fig 4 (a) and (b), shows the natural atomic charges in coulombs (C) for the imidazoline derivative molecules at geometry optimization as well as at 353 K. It is observed that the atoms possess negative charges except for C8 carbon atom in HDM and HDE. This may be due to inductive effect along the carbon chain which is transmitted through the carbon atom and the N4 and N7 nitrogen atoms in HDM and HDE.
It was also observed that the charges become more negative at the structures at 353 K from the geometry optimized structures. The hetero atoms involved (oxygen and nitrogen) are more negative than the carbon atoms except for the last carbon atom at the tail end of the molecules, this may be due to dipole moment (uneven distribution of charges) observed between the last carbon atom and the last hydrogen atom. The oxygen atom was also found to be more negative than the nitrogen atoms. The hydrogen atoms present in the imidazoline molecules are all positively charged.

Bond between atoms in HDM (a)
Geo Opt

Quantum chemistry calculation
Quantum chemical methods are particularly significant in the study of corrosion as it provides researchers with a relatively quick way of studying the structure and behaviour of corrosion inhibitors (Gece 2008, Obot and Obi-Egbedi 2010, Musa et al. 2010, Ju et al. 2008, Arslan et al. 2009, Roque et al. 2008, Masoud et al. 2010, Popova et al. 2007. Quantum chemical parameters such as energy of highest occupied molecular orbital (E HOMO ), energy of lowest unoccupied molecular orbital (E LUMO ),energy gap (ΔE), dipole moment (), energy of deformation (Đ), electron affinity (A), ionization potential (I), global hardness (ƞ), global softness (S)and the Fukui indices were calculated using the Vamp and the Dmol 3 module for the structures at geometry optimization and at the structures at 353 K and the results are shown in Table  4.

HDM HDE
Quantum parameter  Table 4: Calculated quantum chemical parameters for the three molecules studied According to the frontier molecular orbital theory (FMO) of chemical reactivity, transition of electron is due to interaction between highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) of reacting species (Musa et al, 2010a). The electron donating ability of a molecule is usually attributed to E HOMO , which is how well a molecule can give out electrons to the vacant d-orbital of a metal. E LUMO is usually attributed to the ability of a molecule to collect electrons present in a metal's d-orbital, resulting in the creation of another bond (Zarrouk et al, 2012). Therefore, increasing values of E HOMO indicates higher tendency to donate electron(s) to the appropriate acceptor molecule with low energy and empty molecular orbital. Therefore, high values of E HOMO promote adsorption, leading to better inhibition efficiency. From Table 4 HDM has the highest value of E HOMO at geometry optimized molecular structure and at 353 K meaning that these molecules can give out electrons to the vacant d-orbital of a metal more efficiently. The negative signs observed on the values of E HOMO shows that the adsorption is physisorption (Yurt et al, 2006). This is in line with the assumption made concerning the distance observed between the molecules and the Fe surface ( Table 2) that the adsorption observed by the molecules on the Fe surface may be a physical adsorption (physisorption). Lower E LUMO value makes the molecule to easily acquire electrons from the dorbital of a metal surface resulting in the creation of another bond and also to the greater amount of energy released due to the addition of an electron to the molecule (electron affinity). Also from Table 4 HDE has

K
the lowest E LUMO value at geometry optimized molecular structure and at 353 K. The value of ΔE (E LUMO -E HOMO ), the difference in orbital energy between E LUMO as well as E HOMO is a vital stability sign. Higher ΔE values, give good stability and low reactivity and hence bad inhibition efficiency. Low value of ΔE gives better inhibition efficiency, for the reason that lower excitation energy will be required to take away an electron from the last occupied orbital (Khaled et al, 2010). It is seen from Table 4, that HDM has the lowest energy gap at geometry optimized structure and at 353 K. This is in line with the binding energy results obtained at 353 K from Table 2. Fig 5 gives the HOMO and LUMO orbital plot; it can be observed that the HOMO location in the molecules is mostly distributed at the imidazoline ring and few of the plots are observed in the nitrogen and oxygen atoms forming the amine and hydroxyl group joined to the imidazoline ring. This is synonymous to the preferred sites for electrophilic attack (f -) on the molecules. We can say that the part of the molecules with high HOMO density will be oriented towards the mild steel surface and the adsorption could be done by sharing the lone pair of electrons of nitrogen atom on the ring and the  electrons of the aromatic ring with the Fe surface (Hasanov et al, 2007). The LUMO which is synonymous to the preferred sites for nucleophilic attack (f + ) is mainly located at the double bonded carbon atoms inthe alkyl hydrophobic tail in HDM, while HDELUMO orbital is located mainly in the imidazoline ring and the double bonded carbon atom present in the alkyl hydrophobic tail. Note that the blue and yellow isosurfaces depict the electron density difference; the blue regions show electron accumulation, while the yellow regions show electron loss.
The most widely used quantity to describe polarity is the dipole moment of the molecule. Dipole moment is a measure of polarity of a polar covalent bond. It is defined as the product of charge on the atoms and the distance between the two bonded atoms. The dipole moments as well as energy of deformability are parameters characterizing the interaction between molecules (Xia et al, 2008). Molecules with high dipole moment tend to interact with other molecules through electronic interactions (i.e. dipole-dipole interactions). Deformation energy is the energy required to change the orientation of a molecule. The energy of the deformability increases with the increase in dipole moment, making the molecule easier to adsorb on a metal surface. It is shown from the calculation that there was no obvious linear correlation between the values of the dipole moment with the trend of neither the energy gap nor the binding energy and interaction energy as expected. Ionization energy is a fundamental descriptor of the chemical reactivity of atoms and molecules. High ionization energy indicates high stability and chemical inertness and small ionization energy indicates high reactivity of the atoms and molecules. HDM has the least value of ionization energy at geometry optimization and at 353 K. This is also in accordance with the binding energy results. Electron affinity of an atom or molecule is defined as 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 (IUPAC, 1997). HDM has the largest amount of energy release at geometry optimization and at 353 K. This is in line with the binding energy and also the energy gap results (Stupnisek-Lisac et al. 1994, Tang et al. 2006.

HOMO
LUMO f Another important property to measure the molecular stability as well as reactivity of a molecule is the absolute hardness and softness. Chemical hardness fundamentally suggests the resistance towards the deformation or polarization of the electron cloud of the atoms, ions or molecules under small perturbation of a chemical reaction. A hard molecule has a large energy gap and a soft molecule has a small energy gap (Obi-Egbedi et al, 2011). In our present study HDM shows the lowest hardness at geometry optimization and also at 353 K, and shows the highest softness at 353 K. The two molecules show equal softness at geometry optimization. Normally, the inhibitor with the least value of global hardness (hence the highest value of global softness) is expected to have the highest inhibition efficiency (Lukovits et al, 2011). For the simplest transfer of electron, adsorption could occur at the part of the molecule where softness(S), which is a local property, has a highest value (Yurt et al, 2006).
To determine the active site of a molecule, three influencing and controlling factors: natural atomic charge, distribution of frontier molecular orbital and Fukui indices have to be considered. Local reactivity is analyzed by means of the condensed Fukui function. Condensed Fukui functions allow us to distinguish each part of the molecule on the basis of its distinct chemical behavior due to the different substituent functional groups. The nucleophilic and electrophilic attack is controlled by the maximum values of f + and f respectively. The calculated Fukui indices for nucleophilic and electrophilic attack for the two selected inhibitors are tabulated in Table 5 ( only C, N, and O are quoted) and their active sites (Fukui indices plots) are shown on the molecules in Fig 5. It can be seen from Table 5 that the larger values off considering both the geometry optimized and the structures at 353 Kare located on the N4 and N7 atoms for both HDM and HDE, which indicates that these two atoms prefer to form a chemical bond by donation of electrons to the metal surface. With N7 having the highest Fukui indices value, making it the best sites for electrophilic attack for the two molecules studied. While the largest values off + are located on the C19, C8 and also N7 atoms in HDM and C8 atomin HDE, which further suggest that these atoms are responsible to form a feedback bond by the acceptance of electron from the metal surface. But the particular atom which is more susceptible to electrophilic and nucleophilic attack in each of the molecules at geometry optimization and at the structures at 353 K are highlighted in Table 5 (Hosseini et al. 2003, Subramanyam et al. 1993, Domenicano and Hargittai 1992, Levine 1991, Murrrell et al. 1985, Gruber and Bus 1989, Fukui 1975

CONCLUSIONS
In conclusion, the molecular dynamic simulation results indicate that the imidazoline derivatives will most likely use the imidazoline ring to effectively adsorb on the surface of iron, with the alkyl hydrophobic tail forming an n shape at both at geometry optimization and at 353 K. These properties possessed by these molecules (imidazoline derivatives) can prevent water from coming in contact with the Fe surface to a very large extent, based on the n shape canopy-like covering on the Fe surface which is created by the alkyl hydrophobic tail. The quantum chemical calculation based on the natural atomic charge, the frontier molecular orbital and the Fukui indices values and plots shows the active sites of the molecules to be the imidazoline ring mainly the N=C-N region on the ring, the nitrogen and oxygen heteroatoms and the double bonded carbon atoms in the hydrophobic tail of the imidazoline derivative molecules. The bond length analysis also shows the effectiveness of the hetero atoms in terms of reactivity of the bonds on which they share with the carbon atom. The quantum chemistry calculation also reveals that the amine group in HDM and the hydroxyl group in HDE which is attached to the imidazoline ring does not results in a significant increase in the HOMO nor the LUMO density which can aid adsorption. HDM has a lower energy gap of 4.434 eV and 3.824 eV, a higher E HOMO of -4.273 eV and -4.152 eV and a higher global softness of 0.45 and 0.52 compared to HDE which have an energy gap of 4.476 eV and 4.084 eV, a E HOMO of -4.349 eV and -4.607 eV and a global softness of 0.45 and 0.49 at geometry optimization and at 353 K. The adsorption ability of the molecule is given as at geometry optimization HDM > HDE and at 353 K HDM > HDE. Theoretically HDM will be preferred over HDE as a corrosion inhibitor. The adsorption ability of the molecule is very much in line with the binding energy at the temperature studied.