New QSPR model for prediction of corrosion inhibition using conceptual density functional theory

The relationship between structure and corrosion inhibition of a series of twenty-eight quinoline and pyridine derivatives has been established through the investigation of quantum descriptors calculated with PBE/6–311 +  + G** method. A quantitative structure–property relationship (QSPR) model was obtained by examining these descriptors using a genetic algorithm approximation method based on a multiple linear regression analysis. The results indicate that the efficiency of corrosion inhibitors is strongly associated with hardness (η), minimal electrostatic potential (ESPmin), and volume (V) descriptors. Furthermore, the validity of the proposed model is corroborated by an adsorption study on an iron surface Fe(110).


Introduction
Corrosion is the destructive attack of a metal by either chemical or electrochemical processes in the environment [1]. To avoid this process, corrosion inhibitors have been used for several decades; the addition of these chemical substances in low concentrations efficiently decreases the corrosion [2][3][4][5]. In this context, several organic inhibitors have displayed good anticorrosion efficiencies, particularly those that include heteroatoms in their structure as well as aromatic rings [6][7][8][9][10].
On the other hand, several studies have shown that anticorrosion efficiency of chemical compounds can be related to properties and indicators, such as molecular orbital energies, frontier orbital energy gap E LUMO -E HOMO , atomic charges, dipole moments, hardness, chemical potential, electrophilicity, and aromaticity [11][12][13][14][15][16][17][18]. Therefore, a quantitative structural-property relationship (QSPR) can be used to predict the corrosion inhibition of compounds that have not been experimentally synthesized or tested, allowing for the minimization of experimentation time, cost of chemicals, and equipment usage. In this approach, the relationship between either the physicochemical or electronic properties of a selected set of organic compounds and their corrosion is determined [19][20][21][22][23][24][25]. In addition, the aromatic-type molecules are better corrosion inhibitors than the aliphatic ones. The main goal of these studies was to yield suitable mathematical models with predictive potential based on a small set of variables [26][27][28][29][30].
Herein, we describe the chemical reactivity of pyridine and quinoline derivatives and develop a mathematical model for the prediction of the corrosion inhibition by using quantum-chemical descriptors obtained with conceptual DFT methods. In addition, the capability of these chemicals to inhibit the corrosion process has been studied by performing an analysis on an iron surface Fe(110).

Computational procedure
The optimization calculations were performed on twentyeight heterocyclic aromatic molecules based on quinoline and pyridine frameworks. We used the conceptual density functional theory with the Perdew, Burke, and Ernzerhof (PBE) [31, 32] exchange correlation functional and 6-311 + + G** orbital basis set for all atoms using the Gaussian 09 program (Rev. C 01) [33]. The PBE/6-311 + + G** set has yielded reliable and consistent results in a wide variety of chemical systems [34]. All molecules were optimized in gas phase. A second derivative test was applied to each calculation in order to corroborate a minimum on the potential energy surface. A number of chemical descriptors of reactivity such as E HOMO , E LUMO , GAP, hardness (η), electrophilicity (ω), potential ionization (I), dipole momentum (μ D ), molecular electrostatic potential (ESP (max) , ESP (min) ), and aromaticity descriptors anisotropy ( bq ANS 1 ) and isotropy ( bq ISO 1 ) were obtained by means of a single-point calculation in a system formed by a test charge positioned at 1.0 Å from the centroid of the pyridinic ring [35,36]. The harmonic oscillator model of aromaticity index (HOMA), para-delocalization index (PDI), aromatic fluctuation index (FLU), para-linear response (PLR), Shannon index (SA), and the structural parameter molecular volume (V) were obtained with the Multiwfn software [37].
Several mathematical polynomial multiple regression models were obtained for descriptors by using QSARINS software based on the strategy developed by Gramatica [38,39].

Theory and equations
The reactivity of the heterocyclic compounds was analyzed by the following quantum chemical parameters based on conceptual DFT calculations: energy of the highest occupied molecular orbital (E HOMO ), energy of the lowest unoccupied molecular orbital (E LUMO ), and dipolar moment (µ D ). These quantities let us calculate the GAP using Eq. (1).
Additionally, we calculated the vertical ionization potential (I) and vertical electronic affinity (A) [40] by applying Eqs. (2) and (3), where E N−1 c is the system energy without one electron (cation), E N o is the ground state energy (neutral), and E N+1 A is the system energy with one electron added (anion).
The chemical potential and hardness [41][42][43] were obtained from Eqs. (4) and (5), respectively. Similarly, we performed calculations of electrophilicity with Eq. (6) to establish the charge transfer process, i.e., the relationship between the maximum electron transfer and the energy variation. From the equations, it is evident that the chemical reactivity ultimately depends on the electron affinity as well as on the ionization potential.
A Also, the aromaticity descriptors give us a way to measure the electronic delocalization in our systems. The paradelocalization (PDI) and para-linear response (PLR) indexes were used to measure aromaticity of six-membered rings; the larger the indexes, the more aromatic are the systems. The aromatic fluctuation index (FLU) is an aromaticity index based on delocalization among adjacent atoms and can be used to study rings with any number of atoms. Low FLU values correspond to higher aromaticity. The Shannon index (SA) measures aromaticity based on electron density at bond critical points (BCP) in the ring. A low SA value indicates more aromaticity. The 0.003 < SA < 0.005 range is chosen as the boundary of aromaticity/antiaromaticity. The HOMA for a n-membered ring is based on the lengths of individual bonds in comparison with a reference database. If HOMA is equal to 1, it means that the length of each bond is identical to the optimal value and, therefore, the ring is fully aromatic. If HOMA equals to 0, it means that the ring is nonaromatic. If HOMA has a significant negative value, then the ring shows an antiaromaticity characteristic. Finally, the Bird index is another geometry-based quantity used to measure aromaticity; a Bird index value close to 100 means significant aromaticity [34].
We also studied the adsorption of six heterocyclic aromatic molecules on the Fe(110) model systems by means of calculations based in the framework of the DFT method with the generalized gradient adjustment (DFT/GGA), using the exchange correlation functional of Perdew, Burke, and Ernzerhof [31,32]. The effect of the atomic cores on the valence electron density was considered by means of the projected augmented plane-wave (PAW) method of Blöchl [44] as implemented by Kresse and Joubert [45] in the Vienna Ab initio Simulation Package (VASP) [46,47]. This representation of the core states allows to obtain results that converge with a cutoff kinetic energy of 400 eV for the lane wave basis set. The Monkhorst-Pack scheme [48] was used to select the special k-points used to carry out the numerical integrations in the reciprocal space. A conjugated gradient algorithm with an energy criterion of 0.001 eV was used for the atomic convergence, which ensures that forces are smaller than 0.03 eV/Å in all cases.
The Fe(110) surface was characterized using a periodic model of three 6 × 6 unit cell slabs, each having 36 atoms for a total of 108 iron atoms. The periodic slabs in the z-direction were separated using a vacuum region of 16 Å to maintain enough separation from each other. The slabs were constructed using the lattice parameter optimized for bulk Fe(110) and reported in a previous work [49].

Results and discussion
The chemical structures used in this study are presented in Table 1. In order to obtain a reliable model of structure-efficiency relationship on corrosion inhibition, we used the experimental efficiencies of inhibition-to-corrosion data of twenty-eight molecules previously reported; these data were verified under similar experimental conditions [14,50]. Table 1 Chemical structure of studied compounds First, the molecules were optimized in gas phase with PBE/6-311 + + G**. In addition, the calculated structural parameters for this type of heterocyclic compounds, regardless of the functional/basis used, are essentially identical to those reported in the literature [14,14]. Considering these results, we show in Table 1 the structural data for all pyridine and quinoline derivatives obtained at the PBE/6-311 + + G** level of theory.
The reactivity and structural descriptors considered to interpret the efficiency of corrosion inhibition are listed in Table 1 for PBE/6-311 + + G**, while reactivity and aromaticity descriptors are listed in Tables 2 and 3.

QSPR model
Multiple linear regression analyses were used to search the optimal QSPR models that identify the essential parameters involved in corrosion inhibition. The models were validated by means of several statistical parameters [28][29][30], including leave-one-out cross-validation [51] Q 2 LOO , squared correlation coefficient (R 2 ), squared correlation adjusted coefficient (R 2 adj ), the QUIK rule (δK) based on the K multivariate correlation index, and the asymptotic Q 2 rule (δQ).
Moreover, we considered applicability domain (AD), which means that all untested compounds corresponding to the AD of the model will be predictable, and the compounds that are not in the AD area are extrapolations based on the model obtained.
First, when applying the genetic algorithm for selecting reliable molecular descriptors, in this context, the best QSPR model was given by: where %IE is the corrosion inhibition efficiency, the is obtained from Eq. 1 in eV, and ESP min and V were obtained from Multiwfn. This model highlights the correlation   between the %IE with electronic and structural properties for the training set of 28 heterocyclic molecules. The model was obtained with a R 2 = 0.9304 and a Q 2 LOO = 0.9018 . The difference between these values is 0.0286, which is the measure of the dependence of the model on the molecules that are part of it. At first, on the model, a 95% confidence interval was used, obtaining the p values for the variables , ESP min , and V of 0.0000 eV, 0.0319 eV, and 0.0001 Å 3 , respectively. These three p values are smaller than 0.05, so it can be concluded that they are statistically significant.
In Figs. 1 and 2, we show six selected structures as well as the HOMO and LUMO molecular orbitals. These structures were selected to validate the reliability of the model, due to the different corrosion efficiency values they have (i.e., low, medium, and high). Thus, in these figures (Figs. 1 and 2), we distinguish that the delocalization of the electronic densities of HOMO and LUMO orbitals shows that the pyridine and quinoline molecules could donate/accept electrons to/ from the metal d-orbitals through aromatic rings, but their reactivity changes by the effect of the functional group as well as the heteroatom.
In Fig. 3, a comparison is made between the experimental value (experimental endpoint) and the predicted model value, from which it can be deduced that the predicted values have a similar fit to the experimental values. The selected compounds are indicated in the same graph, the blue points correspond to pyridine and its derivatives, and the pink points correspond to quinoline and its derivatives. It is observed that the trend in experimental inhibition is the same as with the proposed model. Figure 4 shows the model residuals, which are a proxy for the regression errors. From a visual analysis of the graph, it can be deduced that the variance is constant (homoscedasticity) and that the residuals do not seem to be correlated. Figure 5 shows the models' residuals, which seem to be homoscedastic and not correlated. The model assumptions on the residuals are met, so it is feasible to use the proposed model.
Additionally, a principal component analysis was performed; we concluded that the variables with the largest proportion of explained variance are , ESP min , and V , confirming that the selected model regressors are representative for obtaining the percentage of corrosion inhibition.

Adsorption of the aromatic molecules on the Fe(110) surface
The adsorption of six heterocyclic aromatic molecules on Fe(110) model systems was also studied. To perform the adsorption studies, three pyridines and three quinolines were selected based on their low, medium, or high reactivity that was observed in our previous results (see Figs. 2 and 6). Each molecule, i.e., pyridine, 2-chloropyridine, 3-benzylpyridine, quinoline, 6-chloroquinoline, and 3-aminoquinoline, was optimized previously and subsequently adsorbed    The adsorption energies obtained by periodic calculations based on DFT are − 48.90, − 83.60, and − 107.02 kcal/mol for pyridine, 2-chloropyridine, and 3-benzylpyridine, respectively, and the energies are − 71.45, − 81.40, and − 126.36 kcal/mol for quinoline, 6-chloroquinoline, and 3-aminoquinoline, respectively. It should be noted that these results show the same trend as the reactivity calculations and structural descriptors presented in the previous section. For example, 3-benzylpyridine turns out to be a better corrosion inhibitor than pyridine. Similarly, the efficiency of 3-aminoquinoline is greater than that of quinoline (see Table 4).

Conclusions
We have developed a reliable model of corrosion inhibition efficiency as function of , , and V which is able to predict, with high accuracy, the corrosion-inhibition attributes of nitrogen-based heterocyclic inhibitors. This model probe is a valuable tool to predict the efficiency of aromatic corrosion inhibitors. Furthermore, the adsorption study corroborates the results because it shows the same trend as the reactivity calculations and structural descriptors.