Theoretical evaluation of the corrosion inhibition performance of 1,3-thiazole and its amino derivatives

2.1 Computational details

Density functional theory (DFT) calculations were carried out using the B3LYP (Becke, 1993) functional, together with the standard double-zeta plus polarization 6-311G (d, p) (Harihara and Pople, 1973) basis set implemented in the GAUSSIAN 03 software. Following the standard nomenclature the calculation will be referred to as B3LYP/6-311G^∗∗, and the frequency analysis was performed to ensure that the researched molecules reached their respective ground state. Generally, the phenomenon of electrochemical corrosion occurs in the liquid phase, so it is necessary to consider the effect of solvent in the calculations. Thus, self-consistent reaction field (SCRF) theory, with Tomasi’s polarized continuum model (PCM) was applied to perform the computations in solution. This method models the solvent (water) as a continuum of uniform dielectric constant (ε = 78.54) and defines the cavity where the solute is placed as a uniform of series of interlocking atomic spheres (Lesar and Milošev, 2009). What needs to stress, local reactivity has been analyzed by means of Fukui Indices (Parr and Yang, 1984), and these properties were calculated with DMol³ version 5.0 available in Materials Studio (Ziegler, 1991), using a BLYP (Becke exchange plus Lee–Yang–Parr correlation) functional (Hamprecht et al., 1998) and a double-numeric quality basis set with polarization functions (DND), which size is comparable with Gaussian 6-311G^∗∗ basis set.

A prediction of iron crystal morphology was necessary before the simulation process. Three representative methods, i.e., the Bravais–Friedel–Donnay–Harker (BFDH) method (Donnay and Harker, 1937; Prywer, 2004), the surface free-energy (SFE) model (Bisker-Leib and Doherty, 2001), and the Attachment Energy (AE) method (Hartman and Perdok, 1955), were utilized to predict crystal structures. More detailed descriptions about these models can be found in Yang and Dong (2011).

The MD simulation of the interaction between the 1,3-thiazole and its amino derivatives molecules and the Fe (1 1 0) surface was carried out in a simulation box (2.97 × 2.97 × 5.05 nm) with periodic boundary conditions in order to model a representative part of an interface devoid of any arbitrary boundary effects. Using six layers of iron atoms gives a sufficient depth that the inhibitor molecules will only be involved in non-bond interactions with iron atoms in the layers of the surface, without increasing the calculation time unreasonably. This also ensured that the depth of the surface was greater than the non-bond cutoff used in calculation. Two types of solvent layers were constructed: one consisted of 1000 H₂O molecules and a single dissolved inhibitor molecule, and the other only consisted of one inhibitor molecule.

The simulation was performed under 298 K, NVT ensemble, with a time step of 1 fs and the simulation time of 200 ps. The temperature was controlled by Andersen thermostat (Andersen, 1980). The initial velocity of each molecule accorded with Maxwell–Boltzmann distribution. Based on hypotheses, such as periodic boundary condition and equivalence of time-average and ensemble-average, Newtonian motion equations were solved with Velocity Verlet algorithm. Van der Waals and coulomb interactions were calculated by atom-based method, and the cutoff radius was 0.95 nm. During the process of simulations, all the atoms in the Fe (1 1 0) surface were fixed. COMPASS (Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies) (Bazooyar et al., 2012) was used to optimize the structures of all components of the system of interest. It is a powerful force field supporting atomistic simulation of condensed phase materials and enables accurate and simultaneous prediction of structural, conformational, vibrational, and thermophysical properties for a broad range of molecules in isolation and in condensed phase, and under a wide range of conditions of temperature and pressure.

2.2 Definitions of global reactivity descriptors

According to density functional theory, for an N-electron system with total energy E and external potential $(\upsilon (\overrightarrow{r}))$ , electronegativity (χ) known as the negative of chemical potential (μ), has been defined as the first derivative of the E with respect to N at $\upsilon (\overrightarrow{r})$ (Parr et al., 1978).

The natural way to approximate the electronic chemical potential μ and the chemical hardness η in DFT is directly evaluating from the ionization potential (IP) and electron affinity (EA) energies (Boshra et al., 2013):

The global electrophilicity index (φ), which measures the electrophilic power of a molecule, has been given the following simple expression (Parr et al., 1999):

3.1 Geometries and global reactivity

The B3LYP/6-311G^∗∗ fully optimized structures of two isomeric forms of amino-1,3-thiazole and of 1,3- thiazole itself are illustrated in Table 1, where selected geometrical parameters in the gas phase and in aqueous solution, as well as the numbering of the atoms of interest necessary for discussion are given. The equilibrium bond lengths (BL), together with the bond orders (BO) for the bonds in the five member ring and the distance from the amino group to the corresponding ring atom are presented in the table. TA shows ideal planar structure of C_s symmetry. In the amino derivatives, the nitrogen atom N6 has the sp³ hybridization. The calculated bond lengths in the gas phase and aqueous phase do not differ significantly. It can be observed that for TA, the BOs are contained in the 1.2–1.6 range. In consistency with the calculated equilibrium bond lengths, these BO values suggest a relatively strong aromatic character for the five-member ring of thiazole. Different BOs were found for 2-ATA. Two BOs are close to 1.1, implying a single bond between the sulfur and carbon (C2 and C5) atoms, while in the C2–N3–C4–C5 region, the BO ranges from 1.2 up to 1.5. This suggests a considerable delocalization of the electronic cloud in this region, that is, in 2-ATA the sulfur atom has decreased considerably the aromatic character of the ring. But we found little evidence this is happening in the 4-ATA case.

Table 1 Selected geometrical parameters and Mayer bond orders of 1,3-thiazole and the amino-1,3-thiazole isomers in gas and aqueous phases.

Bond	Phasea
		TA		2-ATA		4-ATA
		BL (Å)	BO	BL (Å)	BO	BL (Å)	BO
S1–C2	G	1.74783	1.2518	1.77272	1.1653	1.73541	1.2904
	A	1.74344	1.2753	1.77492	1.1702	1.73327	1.3068
C2–N3	G	1.29630	1.6150	1.29975	1.5350	1.29788	1.5940
	A	1.30014	1.5888	1.30618	1.5020	1.30083	1.5689
N3–C4	G	1.37546	1.2918	1.37820	1.2671	1.37934	1.2254
	A	1.37846	1.2781	1.38117	1.2511	1.38327	1.2172
C4–C5	G	1.36250	1.5095	1.35509	1.5607	1.37406	1.3814
	A	1.36234	1.5180	1.35484	1.5723	1.37568	1.3788
C5–S1	G	1.73154	1.2785	1.75504	1.1974	1.73818	1.2608
	A	1.73159	1.2779	1.75638	1.1894	1.73897	1.2564
C2–N6	G			1.36068	1.0748
	A			1.35352	1.1213
C4–N6	G					1.37139	1.0511
	A					1.36875	1.0733

a G, gas phase (ε = 1.0); A, aqueous phase (ε = 78.5).

The effectiveness of an inhibitor can be related to not only its spatial molecular structure, but also with their molecular electronic structure. According to frontier orbital theory, the reaction of reactants mainly occurred on the HOMO and LUMO, and the formation of a transition state is due to an interaction between the frontier orbitals of the reactants. So, it was important to investigate the distribution of HOMO and LUMO for exploration of inhibition mechanism. Organic substances with a higher energy level of HOMO easily donate electrons from HOMO to an empty orbital of appropriate acceptors and E_LUMO denotes the ability of the molecule to accept electrons. As shown in Fig. 1, all the three inhibitor molecules have similar HOMO and LUMO distributions, which were all distributed through the whole molecule. This kind of distribution favored the preferential adsorption of inhibitor on metal surface by following two fashions: one was that the unoccupied d orbitals of Fe atom accepted electrons of inhibitor molecule to form coordinate bond, the other was that the inhibitor molecule accepted electrons from Fe atom with its anti-bonding orbitals to form back-donating bond. Fig. 2 shows the correlation diagram of frontier molecular orbitals for the investigated inhibitors and their calculated HOMO–LUMO gap (ΔE). It was generally acknowledged that low values of ΔE will provide good inhibition efficiency, because the energy for removing an electron from the last occupied orbital will be low (Gece, 2008). In Fig. 2, the values of E_HOMO of three amino acid compounds increased in the following order: 4-ATA > 2-ATA > TA. The values of ΔE decreased in the following order: TA > 2-ATA > 4-ATA. The trend is even more marked in aqueous phase. Thus, 4-ATA has the highest E_HOMO and the lowest ΔE among these three thiazole compounds. Therefore, 4-ATA has the strongest interaction with iron surface and the best inhibition effect in theory.

Close

Figure 1

Schematic representation of HOMO and LUMO molecular orbital of the studied molecules (Color code: N, blue; S, yellow; C, gray; H, white; the same hereinafter).

Export to PPT

Close

Figure 2

Correlation diagram of the frontier molecular orbitals of the inhibitors investigated, and their calculated ΔE (eV) G, gas phase (ε = 1.0); A, aqueous phase (ε = 78.5).

Export to PPT

The effectiveness of 1,3-thiazole and its amino isomers as inhibitor has been further addressed by evaluating the global reactivity parameters. These include global softness (S), electronegativity (χ), global chemical hardness (η), global electrophilicity index (φ), ΔN, and Δψ. The corresponding results are shown graphically in Fig. 3. According to HSAB principle, bulk metals are chemically the softest materials , this would imply that the higher is the S of inhibitor molecules the stronger is their interaction with metal surfaces. On this basis the results of Fig. 2 suggest the following trend of the inhibition effectiveness: 4-ATA > 2-ATA > TA. Furthermore, 4-ATA has the largest fraction of transferred electron to the iron metal, closely followed by 2-ATA, while the TA have the smallest value. The trend in molecule–metal interaction strength, |Δψ|, is also 4-ATA > 2-ATA > TA. There also involves a less important descriptor to these electron-rich inhibitors, i.e., global electrophilicity index; It is generally accepted that the higher the value of φ, the higher the capacity of the molecule to accept electrons. Compared with TA, there exists azyl (–NH₂) grant electron effect in 4-ATA/2-ATA, resulting in a lower value of φ. But when comparing the two isomers, 4-ATA displays greater efficiency in attracting electron ability, and hence the feedback bond will be strengthened. Collectively, the calculated global reactivity parameters follow the order 4-ATA > 2-ATA > TA, in agreement with the above ordering supported by electronic parameters.

Close

Figure 3

Calculated global reactivity parameters of 1,3-thiazole and the amino-1,3-thiazole isomers in gas phase and aqueous phase.

Export to PPT

Based on the data in Table 2, the PCM calculations show that the relative energies of the compounds are lower when increasing the polarity of the solvent. The solvation of the inhibitor molecules causes significant changes in the charge distribution (see Table 3, Mulliken atomic charges q), giving rise to a larger dipole moment than obtained in the gas-phase calculation. In some ways, though, the solvent effect produced small changes in the electronic structure of the inhibitor molecules, but the two phases follow the same trend for the presented properties of three molecules as discussed above. There seems to be no agreement in the literature concerning the correlation between the dipole moment and the inhibition effectiveness (Kokalj, 2010), so we did not employ the molecular dipoles to be used in a blind fashion to judge the inhibition effectiveness.

To describe the energy change associated with the mentioned two processes (i.e., a molecule that is going to receive/back-donate a certain amount of charge), the second order simple charge transfer formula was regarded as a two-parameter expression, in which the donation and back-donation processes are differentiated through the use of the values of the chemical potential for each case, while the hardness is fixed to the value of η = (δ⁺ − δ⁻) in both situations. Thus, according to the simple charge transfer model for donation and back-donation of charges proposed recently by Gomez et al. (2006)) and implemented by (Obi-Egbedi et al., 2011), “when the molecule receives a certain amount of charge, ΔN⁺,

3.2 Local reactivity

The local reactivity of the inhibitors was analyzed through an evaluation of the Fukui indices, which is indicative of the reactive regions and the nucleophilic and electrophilic behavior of the inhibitor molecules. The regions of a molecule where the Fukui function is large are chemically softer than the regions where the Fukui function is small, and 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. For nucleophilic and electrophilic attacks, the condensed Fukui functions, $f_{\text{i}}^{+}$ and $f_{\text{i}}^{-}$ , based on the finite difference approximation (Yang and Mortier, 1986), are given by:

Although atomic charges do not have absolute meaning, they provide a crude but useful description of the electron-density distribution in molecules. In Table 3, the atomic charges derived from Mulliken analysis are given. All the nitrogen and sulfur atoms for each molecule possess an excess of negative charge. The sulfur and nitrogen atoms are assumed to play the key role for the inhibition action because these are the most negative centers of the compounds. Through comparing those Fukui indices in Table 3, it was found that the $f_{\text{k}}^{-}$ at N6 and S1 atoms were the largest for all three molecules, which implied that these atoms preferred to provide electrons to form coordinate bonds with the metal atoms, so these atoms were the electrophilic reactive sites during absorption. As for $f_{\text{k}}^{+}$ , they had the larger value at C2 and C5, which implied that these atoms favored acceptance of electrons from metal atoms to form back-donating. Interestingly, S1 atom has dual characters, i.e., can both donate and back-donate electrons. We hold that this has diminished its adsorption capacity since the donating process is a more important factor in evaluating local reactivity for corrosion inhibitors. Thus, from the Fukui index distribution, it could be concluded that the absorption would occur on multiple reactivity sites, which increased adsorption stability and led to the improvement of inhibition efficiency.

3.3 Selection of adsorption surface

We calculated the surface energy (E_surf) in a slab model as (Arya and Carter, 2004):

Close

Figure 4

The crystal habits of iron simulated by three different models.

Export to PPT

Just for an added delineation, the XRD pattern of α-Fe was obtained by using the Powder Diffraction function in Reflex module included in software Materials Studio, and shown in Fig. 5. In general, the half width (FWHM) of X-ray diffraction is used as a strength index of material surface. It can be seen from Fig. 5 that the value of FWHM increases gradually as diffraction angle increases; the (1 1 0) diffraction peak possesses good intensity and then high orientation and crystallinity.

Close

Figure 5

Simulated XRD pattern of α-Fe (Inset: FWHM vs 2theta) (X-ray, copper radiation, step size: 0.05°).

Export to PPT

In brief, the density packed (1 1 0) surface has the most stabilization, so we choose Fe (1 1 0) surface to simulate the adsorption process.

3.4 Molecular dynamics simulations

The molecular dynamics simulations were performed to study the adsorption behavior of the TA, 2-ATA, and 4-ATA on the Fe (1 1 0) surface. Whether the model system has reached equilibrium or not was ascertained by the equilibrium criterions of temperature and energy simultaneously (Zeng et al., 2011), i.e., the fluctuations of temperature and energy should be confined to 5–10%. Taking 4-ATA as an example, according to the temperature and energy curves from MD simulations at 298 K (Fig. 6), the temperature fluctuates in a range of (298 ± 20) K and the fluctuation of energy is less than 0.5%, indicating that the system has reached an equilibrium state.

Close

Figure 6

Temperature equilibrium curve (a) and energy fluctuant curve (b) of 4-ATA molecule on Fe (1 1 0) surface (results in vacuo).

Export to PPT

The equilibrium configurations for the three corrosion inhibitors are shown in Fig. 7. By careful examination of Fig. 7, it can be noticed that, during the process of simulation, 1,3-thiazole and its derivatives molecules moved gradually near the Fe (1 1 0) surface. After the system reached equilibrium, all the investigated inhibitors adsorbed nearly parallel to the Fe (1 1 0) surface, and they adsorb on the metal surface through the nitrogen (N) and sulfur (S) atoms as well as the five-membered aromatic rings. Quantitative appraisal of the interaction between each molecule and the Fe surface was achieved by calculating the adsorption energy (E_ads) using Eq. (14) (Znini et al., 2012)

Close

Figure 7

Representative snapshots of TA, 2-ATA, and 4-ATA on the metal surfaces, viewed along different perspectives (inset, on-top view) to support a spatial impression.

Export to PPT

To get proper visualization of the surface orientation of the inhibitor molecules, we have modified the Fe slabs by reducing the diameters of the sub-surface Fe atoms (see Fig. 7 insets). In this case, it is especially obvious that those polarizable atoms (C, N, S) along the molecular backbone seem to avoid contact with the Fe atoms on the surface plane {larger spheres on the Fe (1 1 0) slab} and appear to align with vacant sites on the fcc lattice atop the metal surface. In another word, the molecular backbone of the adsorbate molecule is preferentially accommodated in characteristic epitaxial grooves on the metal surface. Such epitaxial adsorption orientations are associated with a minimum free energy of adsorption as the adsorbates can be considered as extensions of the fcc lattice and adsorption strength scales with improved fit of the polarizable atoms of a molecule to multiple epitaxial sites (Heinz et al., 2009). This phenomenon accounts for the high adsorption strengths of the adsorbent/adsorbate systems and is commonly seen in other corrosion inhibitors (Oguzie et al., 2011; Mejeha et al., 2012).

In order to simulate the realistic situation, it is important to take the effect of solvent into account. The corresponding simulation result is presented in Fig. 8, and the calculated interaction energy and binding energy are also listed in Table 5. It can be seen from Fig. 8 that the optimization configurations of 1,3-thiazole and its amino derivatives in aqueous solution are similar to those in the vacuum systems. Moreover, the order of E_ads obtained from calculations in the aqueous solution is the same as the one in the nonaqueous system. Nevertheless, the values of E_binding in aqueous solution are slightly higher than those in the nonaqueous system. The parallel adsorption configurations under different environments ensure that the iron surface can be maximally covered by the investigated inhibitor molecules and form a barrier layer between the iron surface and the aggressive media, thus enhancing the inhibition efficiency.

Close

Figure 8

Equilibrium adsorption configurations of the three inhibitors on Fe (1 1 0) surface in water solution (Inset, on-top view).

Export to PPT

After analyzing the MD simulation results of corrosion inhibitors with iron crystal surface, pair correlation functions of the super-molecular structures, g(r), were obtained (see Fig. 9). As a basic rule, the peak in the g(r) ∼ r curve within 0.35 nm is caused by chemical bonds, and that outside 0.35 nm is caused by coulomb and Van der Waals interactions (Zeng et al., 2011). As shown in Fig. 9, the distance between the 4-ATA/2-ATA and the iron atoms of the Fe (1 1 0) face is 0.325 nm, less than 0.350 nm, which indicates that chemical bonds have formed between these two corrosion inhibitors and iron atoms. However, the Fe-TA distance (0.355) slightly exceeds 0.350 nm, suggesting a weak chemisorption. Furthermore, the countless decentralized peaks located outside 0.35 nm, which are derived from physical interactions, may also contribute to the net molecule surface attraction.

Close

Figure 9

Pair correlation functions of the three inhibitors on Fe (1 1 0) surface in water solution.

Export to PPT