Furfural Analogs as Sustainable Corrosion Inhibitors—Predictive E ﬃ ciency Using DFT and Monte Carlo Simulations on the Cu(111), Fe(110), Al(111) and Sn(111) Surfaces in Acid Media

: Nowadays, theoretical calculation tools have become powerful in predicting the behavior of corrosion inhibitors on the surface of metals and, therefore, avoiding energy consumption and the cost of experimental tests. This work aims to predict the inhibitory power of some furan derivatives on Cu (111), Fe (110), Al (111) and Sn (111) surfaces in acidic media. For this purpose, three furan derivatives—furan-2-carbaldehyde (FF1), 5-(hydroxymethyl)furfural (FF2) and 5-(hydroxymethyl)furoic acid (FF3)—have been selected to compare their intrinsic properties against corrosion as well as their behavior on iron (Fe), copper (Cu), aluminum (Al) and tin (Sn) surfaces in acid medium. Typically, the anti-corrosive properties of FF1, FF2 and FF3 were studied by using quantum chemical calculations and Monte Carlo simulations. Density Functional Theory (DFT), lowest unoccupied (E LUMO ) and highest occupied (E HOMO ) molecular orbital energies, energy gap ( ∆ E), chemical hardness ( η ), softness ( σ ), electronegativity ( χ ), electrophilicity ( ω ) and nucleophilicity ( ε ) have been calculated and discussed. Theoretical vibrational spectra were also calculated to exhibit the functional groups in the selected chemicals. On the other hand, the adsorption behaviors of FF1, FF2 and FF3 were studied on the Fe(110), Cu(111), Al(111) and Sn(111) surfaces. As a result, the adsorption energies of all molecules are ordered as Fe(110) < Cu(111) < Al(111) < Sn(111) and FF3 seems to be more e ﬀ ective as a corrosion inhibitor due to the existence of both carboxylic acid and hydroxyl groups, which consist of favorable sites of adsorption into the metal surface.


Introduction
Corrosion still has a huge economic impact in most industrial countries today, accounting for about 3-4 percentage points of the gross domestic product (GDP) [1][2][3]. As fossil resources are being depleted at a rapid rate, there is an urgent need to develop renewable energies to replace fossil fuels derived from oil. Biomass, the only renewable resource of organic carbon in nature, is a response to oil substitution [1,2]. However, the management of corrosion in biorefinery complexes has only recently been identified as a potential challenge in the scientific literature [4]. From an industrial point of view, the effects related to corrosion must be anticipated in order to apprehend an appropriate evaluation of the possibilities of Capital Expenditure (CAPEX) and Operational Expenditure (OPEX), return on investment and extension of the life of the installations. Indeed, the aging, maintenance and modernization operations of biorefineries can induce new potential corrosive environments that must be identified and reassessed regularly. Operational disruptions, extra costs, industrial risks and environmental impacts can be effective due to the erroneous or underestimated evaluation of corrosion issues. Several studies have evaluated the annual direct cost of corrosion to be about 3.1% of gross national product (GNP) in industrialized countries [1][2][3].
Currently, computational software is becoming a much-used and trusted tool to explain the behavior of corrosion inhibitors in several media and metal surfaces [6,[33][34][35][36]. Density functional theory (DFT) is a computational modeling method usually used to investigate the intrinsic properties of molecules. For inhibitors, it is mainly conducted to predict chemical properties such as the highest occupied molecular orbital-lowest unoccupied molecular orbital (HOMO-LUMO) energy gap, chemical hardness, softness, electronegativity, chemical potential, proton affinity, electrophilicity and nucleophilicity of chemical species. On the other hand, Monte Carlo simulations were used to understand the behavior of inhibitors on the metal surfaces. Indeed, adsorption energy and binding energy are calculated and investigated [37][38][39][40][41][42][43][44][45]. To the best of our knowledge, furanic derivatives were studied using theoretical tools only for conversion [46,47], electrochemical oxidation [48] and hydrodeoxygenation [49,50]. In our work, theoretical studies were investigated to predict the efficiency of three selected furan derivatives: furan-2-carbaldehyde (FF1), 5-(hydroxymethyl) furfural (FF2) and 5-(hydroxymethyl)furoic acid (FF3). Firstly, the intrinsic properties such as the E HOMO -E LUMO energy gap (∆E), chemical hardness (η), electronegativity (χ), the fraction of electrons transferred (∆N), total negative charges and dipole moment of FF1, FF2 and FF3 were studied by using quantum chemical calculations. On the other hand, the adsorption and binding energies of these three molecules were examined on the surface of several metals, namely Cu(111), Fe(110), Al(111) and Sn(111) using Monte Carlo simulations.

Monte Carlo Simulations
In the recent years, Monte Carlo simulations have become a potent computation tool that is employed widely to examine the possible inhibitor-metal interactions for a large number of inhibition systems [41,57]. For this purpose, the interaction between the investigated furfural compounds and selected metal surfaces (Cu, Fe, Al and Sn) was studied using Monte Carlo simulations, as the temperature of the system was slowly decreased (i.e., simulated annealing algorithm [58,59]). The organic inhibitors are wieldy used in the acidic environment, which is considered in the current study. In such an environment, the surface of the selected substrates consists of pure metal [60]. To predict the most adequate (h k l) Miller index of each metal surface for constructing of simulation box, the Bravais-Friedel-Donnay-Harker (BFDH) method was utilized [53]. Five layers and a vacuum region of 60 Å were used to model the slab of each studied metal in the current work. The Van der Waals and electrostatic interactions were calculated via the atom-based and Ewald summation methods, respectively. Materials Studio 6.0 software was utilized to perform these calculations with the Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) force field. These simulations were carried out using five heating ramps with 50,000 steps each. The obtained conformers were pre-geometrized and inputted with a temperature range between 100 and 10 5 K. To get results with excellent quality, the convergence threshold was fixed at 5 × 10 −5 Å, 5 × 10 −3 kcal mol −1 Å −1 and 10 −4 kcal mol −1 for displacement, force and energy, respectively. The charges used in the calculation were given by the employed force field. The adsorption energy (E ads ) of the furfural derivative on each metal surface type was calculated according to the following expression [8]: where E T denotes the total energy of the whole system, E Surf is the energy of the metal surface and E Furfural ads is the energy of the adsorbed furfural derivative on metal surface. In the present work, we have limited the simulations process in the gas phase to reduce the required time for calculations. Additionally, it was reported [61,62] that the aqueous phase can influence the magnitude of calculated adsorption energies of inhibitor molecules, while the trend of these energies remains the same one either with or without solvent.

Optimized Structures
The furfural derivatives FF1, FF2 and FF3 contain several functional groups such as formyl, hydroxyl and carboxylic groups ( Table 1). The optimized structures of inhibitor molecules in aqueous Sustainability 2020, 12, 3304 4 of 14 and gas phases using the hybrid DFT functional method (i.e., B3LYP/6-31G++(2d,p)) are shown in Figures 1 and 2, respectively.  Before launching the other theoretical calculations, the structure optimization step of the investigated molecules was done in aqueous and gas phases to determine the geometrical structure with a local minimum on the conformational. Namely, the precise theoretical determination of the geometry of FF1, FF2 and FF3 chemicals were studied at the minima of their potential energy [33]. Thus, it is essential to obtain optimized structures of studied inhibitors to subsequently determine the other parameters, such as a map of electrostatic potential (ESP).

Electrostatic Potential (ESP) 2D Maps
By using ESP maps, it is possible to visualize the electron distribution, and therefore identify the centers or the areas of their concentration in each compound [63]. The contour maps of electron density reveal that oxygen atoms on the studied molecule inhibitors exhibit favorable interaction sites, taking into account the difference between the oxygen atoms in the functions and the other atom. The oxygen atoms, in this case, exist in the formyl group of FF1, in the formyl and hydroxyl groups of FF2 and in the hydroxyl and carboxylic groups in FF3. Interaction sites surrounded by a dark red contour contribute to form the bonding interactions between metal surfaces and inhibitors [13]. The dark red color in the contour map of negative potential particularly surrounds oxygen molecules and its regions in FF1, FF2 and FF3 molecules in gas and aqueous phases, whereas the green color is scattered in the positive potential region according to the contour representation of electrostatic at the DFT/B3LYP/6-31++G(2d,p) calculation level ( Figure 4). Similar geometries were obtained in the aqueous phase ( Figure 5).   Before launching the other theoretical calculations, the structure optimization step of the investigated molecules was done in aqueous and gas phases to determine the geometrical structure with a local minimum on the conformational. Namely, the precise theoretical determination of the geometry of FF1, FF2 and FF3 chemicals were studied at the minima of their potential energy [33]. Thus, it is essential to obtain optimized structures of studied inhibitors to subsequently determine the other parameters, such as a map of electrostatic potential (ESP).

Electrostatic Potential (ESP) 2D Maps
By using ESP maps, it is possible to visualize the electron distribution, and therefore identify the centers or the areas of their concentration in each compound [63]. The contour maps of electron density reveal that oxygen atoms on the studied molecule inhibitors exhibit favorable interaction sites, taking into account the difference between the oxygen atoms in the functions and the other atom. The oxygen atoms, in this case, exist in the formyl group of FF1, in the formyl and hydroxyl groups of FF2 and in the hydroxyl and carboxylic groups in FF3. Interaction sites surrounded by a dark red contour contribute to form the bonding interactions between metal surfaces and inhibitors [13]. The dark red color in the contour map of negative potential particularly surrounds oxygen molecules and its regions in FF1, FF2 and FF3 molecules in gas and aqueous phases, whereas the green color is scattered in the positive potential region according to the contour representation of electrostatic at the DFT/B3LYP/6-31++G(2d,p) calculation level ( Figure 4). Similar geometries were obtained in the aqueous phase ( Figure 5).   Before launching the other theoretical calculations, the structure optimization step of the investigated molecules was done in aqueous and gas phases to determine the geometrical structure with a local minimum on the conformational. Namely, the precise theoretical determination of the geometry of FF1, FF2 and FF3 chemicals were studied at the minima of their potential energy [33]. Thus, it is essential to obtain optimized structures of studied inhibitors to subsequently determine the other parameters, such as a map of electrostatic potential (ESP).

Electrostatic Potential (ESP) 2D Maps
By using ESP maps, it is possible to visualize the electron distribution, and therefore identify the centers or the areas of their concentration in each compound [63]. The contour maps of electron density reveal that oxygen atoms on the studied molecule inhibitors exhibit favorable interaction sites, taking into account the difference between the oxygen atoms in the functions and the other atom. The oxygen atoms, in this case, exist in the formyl group of FF1, in the formyl and hydroxyl groups of FF2 and in the hydroxyl and carboxylic groups in FF3. Interaction sites surrounded by a dark red contour contribute to form the bonding interactions between metal surfaces and inhibitors [13]. The dark red color in the contour map of negative potential particularly surrounds oxygen molecules and its regions in FF1, FF2 and FF3 molecules in gas and aqueous phases, whereas the green color is scattered in the positive potential region according to the contour representation of electrostatic at the DFT/B3LYP/6-31++G(2d,p) calculation level ( Figure 4). Similar geometries were obtained in the aqueous phase ( Figure 5).  aqueous and gas phases using the hybrid DFT functional method (i.e., B3LYP/6-31G++(2d,p)) are shown in Figures 1 and 2, respectively. In order to highlight the functional groups in FF1, FF2 and FF3, the vibrational spectroscopy of all molecules was calculated using the B3LYP/6-31++G(2d,p) method ( Figure 3). The vibrational spectra of FF1, FF2 and FF3 confirmed the existence of aldehyde, aldehyde and alcohol and carboxylic acid and alcohol, respectively. The peaks at around 3000 cm −1 are mainly attributed to hydroxyl groups. Moreover, the peaks at 1705, 1698 and 1745 cm −1 are attributed to carbonyl groups.  aqueous and gas phases using the hybrid DFT functional method (i.e., B3LYP/6-31G++(2d,p)) are shown in Figures 1 and 2, respectively. In order to highlight the functional groups in FF1, FF2 and FF3, the vibrational spectroscopy of all molecules was calculated using the B3LYP/6-31++G(2d,p) method ( Figure 3). The vibrational spectra of FF1, FF2 and FF3 confirmed the existence of aldehyde, aldehyde and alcohol and carboxylic acid and alcohol, respectively. The peaks at around 3000 cm −1 are mainly attributed to hydroxyl groups. Moreover, the peaks at 1705, 1698 and 1745 cm −1 are attributed to carbonyl groups.  In order to highlight the functional groups in FF1, FF2 and FF3, the vibrational spectroscopy of all molecules was calculated using the B3LYP/6-31++G(2d,p) method ( Figure 3). The vibrational spectra of FF1, FF2 and FF3 confirmed the existence of aldehyde, aldehyde and alcohol and carboxylic acid and alcohol, respectively. The peaks at around 3000 cm −1 are mainly attributed to hydroxyl groups. Moreover, the peaks at 1705, 1698 and 1745 cm −1 are attributed to carbonyl groups.
Before launching the other theoretical calculations, the structure optimization step of the investigated molecules was done in aqueous and gas phases to determine the geometrical structure with a local minimum on the conformational. Namely, the precise theoretical determination of the geometry of FF1, FF2 and FF3 chemicals were studied at the minima of their potential energy [33]. Thus, it is essential to obtain optimized structures of studied inhibitors to subsequently determine the other parameters, such as a map of electrostatic potential (ESP).

Electrostatic Potential (ESP) 2D Maps
By using ESP maps, it is possible to visualize the electron distribution, and therefore identify the centers or the areas of their concentration in each compound [63]. The contour maps of electron density reveal that oxygen atoms on the studied molecule inhibitors exhibit favorable interaction sites, taking into account the difference between the oxygen atoms in the functions and the other atom. The oxygen atoms, in this case, exist in the formyl group of FF1, in the formyl and hydroxyl groups of FF2 and in the hydroxyl and carboxylic groups in FF3. Interaction sites surrounded by a dark red contour contribute to form the bonding interactions between metal surfaces and inhibitors [13]. The dark red color in the contour map of negative potential particularly surrounds oxygen molecules and its regions in FF1, FF2 and FF3 molecules in gas and aqueous phases, whereas the green color is scattered in the positive potential region according to the contour representation of electrostatic at the DFT/B3LYP/6-31++G(2d,p) calculation level (Figure 4). Similar geometries were obtained in the aqueous phase ( Figure 5).  Before launching the other theoretical calculations, the structure optimization step of the investigated molecules was done in aqueous and gas phases to determine the geometrical structure with a local minimum on the conformational. Namely, the precise theoretical determination of the geometry of FF1, FF2 and FF3 chemicals were studied at the minima of their potential energy [33]. Thus, it is essential to obtain optimized structures of studied inhibitors to subsequently determine the other parameters, such as a map of electrostatic potential (ESP).

Electrostatic Potential (ESP) 2D Maps
By using ESP maps, it is possible to visualize the electron distribution, and therefore identify the centers or the areas of their concentration in each compound [63]. The contour maps of electron density reveal that oxygen atoms on the studied molecule inhibitors exhibit favorable interaction sites, taking into account the difference between the oxygen atoms in the functions and the other atom. The oxygen atoms, in this case, exist in the formyl group of FF1, in the formyl and hydroxyl groups of FF2 and in the hydroxyl and carboxylic groups in FF3. Interaction sites surrounded by a dark red contour contribute to form the bonding interactions between metal surfaces and inhibitors [13]. The dark red color in the contour map of negative potential particularly surrounds oxygen molecules and its regions in FF1, FF2 and FF3 molecules in gas and aqueous phases, whereas the green color is scattered in the positive potential region according to the contour representation of electrostatic at the DFT/B3LYP/6-31++G(2d,p) calculation level (Figure 4). Similar geometries were obtained in the aqueous phase ( Figure 5).

HOMO and LUMO Energies and Derived Parameters
The highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) of a chemical molecule are very important in defining its reactivity [64,65]. Figures 6 and 7

HOMO and LUMO Energies and Derived Parameters
The highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) of a chemical molecule are very important in defining its reactivity [64,65]. Figures 6 and 7 represent the HOMO and LUMO orbitals of the studied compounds in gas and aqueous phases, respectively. The compounds FF1, FF2 and FF3 present a clear contribution of ρ-orbitals at the cyclical level (furan ring), consisting of an aromatic ring with five atoms, including one oxygen atom. Figure 5. The contour representation of electrostatic potential regions of negative (positive) potential is red (green) for inhibitor molecules in aqueous phase using the DFT/B3LYP/6-31++G(2d,p) method.

HOMO and LUMO Energies and Derived Parameters
The highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) of a chemical molecule are very important in defining its reactivity [64,65]. Figures 6 and 7 represent the HOMO and LUMO orbitals of the studied compounds in gas and aqueous phases, respectively. The compounds FF1, FF2 and FF3 present a clear contribution of ρ-orbitals at the cyclical level (furan ring), consisting of an aromatic ring with five atoms, including one oxygen atom. Figure 6. The highest occupied molecular orbitals (HOMOs) and lowest unoccupied molecular orbitals (LUMOs) of inhibitor molecules in gas phase using the DFT/B3LYP/6-31++G(2d,p) method. Figure 6. The highest occupied molecular orbitals (HOMOs) and lowest unoccupied molecular orbitals (LUMOs) of inhibitor molecules in gas phase using the DFT/B3LYP/6-31++G(2d,p) method. Figure 5. The contour representation of electrostatic potential regions of negative (positive) potential is red (green) for inhibitor molecules in aqueous phase using the DFT/B3LYP/6-31++G(2d,p) method.

HOMO and LUMO Energies and Derived Parameters
The highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) of a chemical molecule are very important in defining its reactivity [64,65]. Figures 6 and 7 represent the HOMO and LUMO orbitals of the studied compounds in gas and aqueous phases, respectively. The compounds FF1, FF2 and FF3 present a clear contribution of ρ-orbitals at the cyclical level (furan ring), consisting of an aromatic ring with five atoms, including one oxygen atom. Figure 6. The highest occupied molecular orbitals (HOMOs) and lowest unoccupied molecular orbitals (LUMOs) of inhibitor molecules in gas phase using the DFT/B3LYP/6-31++G(2d,p) method. The quantum molecular results summarized in Tables 2 and 3 aims to describe in detail the energetic and structural characteristics of the studied molecules. The quantum chemical parameters calculated using the HF/6-31G++(2d,p), MP2/6-31G++(2d,p) and B3LYP/6-31G++(2d,p) methods for the inhibitors in aqueous and gas phases are E LUMO , E HOMO , ∆E, I, A, η, χ, ω, ε and ∆N. The E LUMO and E HOMO values are slightly different with a negligible difference for each method, which means that these are representative results. These two elements, E LUMO and E HOMO , are the basis for calculating all the other parameters that are estimated following the Koopmans theorem [66] to determine the properties related to the reactivity and selectivity of the inhibitor. Starting with the energy gap ∆E (i.e., E HOMO -E LUMO ) is an essential parameter as the reactivity of the inhibitor compound, and thus an efficient inhibitor is characterized by a small energy gap [12]. All values of ∆E determined for the inhibitors studied by the different methods have the same order FF1 > FF2 > FF3, which shows that FF3 has good inhibitory power compared to the other inhibitors FF2 and FF1, respectively (Tables 1 Sustainability 2020, 12, 3304 7 of 14 and 2). Besides, the reactivity of the inhibitors can be also estimated by the other parameters, such as electronegativity expressed from E LUMO and E HOMO using Equation 2, depending on the results that show that FF3 has a small electronegativity value compared to FF2 and FF1. Therefore, concerning the reactivity of these compounds, FF3 tends more to react as a donor of electrons [14]. The global hardness (η) values were calculated using Equation 3, with softness values defined using Equation 3; after that, the global electrophilicity index (ω) were expressed in terms of the electronegativity and global hardness parameters according to Equation 3. All of these parameters are essential to deeply describe the reactivity of an inhibitor in addition to the ∆N as determined by Equation 4, which considers the most critical parameter because it combines several previous parameters; thus it gives a global description of the reactivity of an inhibitor. It is remarkable that the global hardness is listed in ascending order as following FF3, FF2 and FF3, and that means that FF3 is more reactive than both FF2 and FF1 (Table 3). This is because it is demonstrated that the chemical hardness of molecules is precisely the screened interaction energy of the electrons in the frontier (HOMO and LUMO) orbitals, which are the resistance of an atom to a charge transfer [67]. For the electrophilicity, which describes a good electrophile, while a small value of electrophilicity describes a good nucleophile that likes positive charges, the latter generally represent the charges of a metallic surface; therefore according to the obtained values, FF3 > FF2 > FF1 in this order react with a Sustainability 2020, 12, 3304 8 of 14 positively charged metal surface. In addition to the number fraction of electrons transferred, ∆N for FF3, FF2 and FF1 are in agreement with all adopted bases, as reported in Table 3, which also confirms that FE3 has a high number of electron transfer, hence predicting that the FF3 has the highest inhibition performance compared to the other inhibitors FF2 and FF1. The number of transferred electrons (∆N) gives information about the number of electrons transferred to the acceptor surface [68].

Charge Distribution
The presence of a high negative charge in some atoms or sites in compounds gives them a higher tendency to donate electrons to react with the metal surface. Knowing that, the reaction sites of the FF1, FF2 and FF3 molecules are the oxygen atoms, namely that of the molecular functions aldehyde ((H-C=O) and alcohol (-OH alcohol) or (-OH acid)). For this reason, one of the most common methods to present this property is Mulliken analysis [16]. The Mulliken partial charges on the different atoms of the optimized molecules studied are shown in Figures 8 and 9 and summarized in Table 4. These Mulliken charge distributions for the molecules are calculated at the B3LYP/6-31G++(2d,p) in the gas and aqueous phase. These distribution charges on some heteroatoms such as oxygen (O) and carbon (C) can make such groups into susceptible active centers, which explains the highest negative charges of some heteroatoms, especially the oxygen atom (O), because it is a part of inhibitors function groups. These groups part, react, or adsorb on a metallic surface.

Charge Distribution
The presence of a high negative charge in some atoms or sites in compounds gives them a higher tendency to donate electrons to react with the metal surface. Knowing that, the reaction sites of the FF1, FF2 and FF3 molecules are the oxygen atoms, namely that of the molecular functions aldehyde ((H-C=O) and alcohol (-OH alcohol) or (-OH acid)). For this reason, one of the most common methods to present this property is Mulliken analysis [16]. The Mulliken partial charges on the different atoms of the optimized molecules studied are shown in Figures 8 and 9 and summarized in Table 4. These Mulliken charge distributions for the molecules are calculated at the B3LYP/6-31G++(2d,p) in the gas and aqueous phase. These distribution charges on some heteroatoms such as oxygen (O) and carbon (C) can make such groups into susceptible active centers, which explains the highest negative charges of some heteroatoms, especially the oxygen atom (O), because it is a part of inhibitors function groups. These groups part, react, or adsorb on a metallic surface.
Furthermore, all negatives atoms on the molecular skeleton can contribute to this process through an intramolecular synergistic effect [15].

Morphology of Studied Metal Surfaces
The obtained BFDH morphology results for selected metal subtracts in the current study were summarized (Table 5). It can be seen that the multiplicities of Cu(111), Fe(110), Al(111) and Sn(111) faces are 8, 12, 8 and 8, respectively. This indicates that these faces have more contact sites to interact with furfural molecules. Furthermore, the larger dhkl distances were obtained in the case of the Cu(111), Fe(110), Al(111) and Sn(111) faces, which reveals that these faces are densely packed metal surfaces. Furthermore, the higher % total facet area values are noted for the Cu(111), Fe(110), Al(111) and Sn(111) faces, which accounts for 78%, 64%, 78% and 100% of the crystal surface, respectively;  Furthermore, all negatives atoms on the molecular skeleton can contribute to this process through an intramolecular synergistic effect [15].

Morphology of Studied Metal Surfaces
The obtained BFDH morphology results for selected metal subtracts in the current study were summarized (Table 5). It can be seen that the multiplicities of Cu(111), Fe(110), Al(111) and Sn(111) faces are 8, 12, 8 and 8, respectively. This indicates that these faces have more contact sites to interact with furfural molecules. Furthermore, the larger d hkl distances were obtained in the case of the Cu(111), Fe(110), Al(111) and Sn(111) faces, which reveals that these faces are densely packed metal surfaces. Furthermore, the higher % total facet area values are noted for the Cu(111), Fe(110), Al(111) and Sn(111) faces, which accounts for 78%, 64%, 78% and 100% of the crystal surface, respectively; this means there are the main faces where the adsorption process of furfural derivatives can occur. Accordingly, the (110), (110), (111) and (111) faces were chosen as representative surface models for copper, iron, aluminum and tin metal, respectively, during the Monte Carlo simulations.

Conclusions
In summary, DFT calculation with different basis sets and Monte Carlo simulations were carried out to compare and evaluate the corrosion inhibition efficiencies of three furan derivatives, which are furan-2-carbaldehyde, 5-(hydroxymethyl) furfural and 5-(hydroxymethyl) furoic acid on Fe(110), Cu(111), Al(111) and Sn(111) surfaces. The theoretical vibrational study confirmed that all of these compounds are distinguished by the functional groups (aldehyde and carboxylic groups). The neutral and portioned forms were optimized and investigated. Using the results of DFT calculations, it was shown that FF3 has a small electronegativity value compared to FF2 and FF1. Therefore, concerning the reactivity of these compounds, FF3 has a larger tendency to react as a donor of electrons. Additionally, Monte Carlo simulations showed that the adsorption energies of all molecules are ordered as Fe(110) < Cu(111) < Al(111) < Sn(111). In conclusion, 5-(hydroxymethyl)furoic acid could be a good candidate for anticorrosive protection and should to be tested and studied, especially on iron surfaces in acidic media.