DFT Study of Azole Corrosion Inhibitors on Cu2O Model of Oxidized Copper Surfaces: I. Molecule–Surface and Cl–Surface Bonding

The adsorption of three simple azole molecules—imidazole, triazole, and tetrazole—and Cl on various sites of several Cu 2 O(111)- and Cu 2 O(110)-type surfaces, including Cu and O vacancies, was characterized using density functional theory (DFT) calculations; the three molecules can be seen as models of azole corrosion inhibitors and Cl as a corrosion activator. Both non-dissociative and dissociative adsorption modes were considered for azole molecules; the latter involves the N–H bond cleavage, hence we also addressed the adsorption of H, which is a co-product of the dissociative adsorption. We find that molecules and Cl bind much stronger to unsaturated Cu sites compared to saturated ones. Dissociated molecules bind considerably stronger to the surface compared to the intact molecules, although even the latter can bind rather strongly to specific unsaturated Cu sites. Bader analysis reveals that binding energies of dissociated molecules at various Cu sites correlate with Bader charges of Cu ions before molecular adsorption, i.e., the smaller the Cu charge, the stronger the molecular bonding. All three azole molecules display similar non-dissociative adsorption energies, but significant difference between them appears for dissociative adsorption mode, i.e., dissociated triazole and tetrazole bind much stronger than dissociated imidazole because the former two can form two strong N–Cu bonds, but imidazole cannot due to its incompatible molecular geometry. Dissociative adsorption is consequently favorable only for triazole and tetrazole, but only at oxygen vacancy sites, where it proceeds barrierlessly (or almost so). This observation may suggest that, for imidazole, only the neutral form, but, for triazole and tetrazole, also their deprotonated forms are the active species for inhibiting corrosion under near neutral pH conditions, where copper surfaces are expected to be oxidized. As for the comparison with the Cl–surface bonding, the calculations indicate that only dissociated triazole and tetrazole bind strong enough to rival the Cl–surface bonds.


Introduction
Some organic molecules are well known for their ability to slow down the corrosion of metals. Among them, azoles-five-membered heterocyclic molecules containing one or more nitrogen atoms-are well known as corrosion inhibitors for copper [1][2][3][4]. Although the mechanism by which organic molecules inhibit corrosion is usually not known, it is nevertheless widely accepted that strong adsorption of inhibitor molecules represents an important step in achieving the inhibitory effect. Indeed, Bockris stated that organic molecules must be adsorbed to become inhibitors [5]. From this point of view, it is therefore important to characterize the molecule-surface bonding in the context of corrosion inhibition. It is here where density functional theory (DFT) based first-principle modeling, utilizing the periodic slab representation of the surface, can provide useful information and consequently helps opening a way towards a more rational design of new corrosion inhibitors [6][7][8][9]. It should be noted, though, that the inhibitor-surface bonding itself is far from synonymous with inhibition of corrosion. Instead, it represents only one aspect towards an atomic-scale understanding of corrosion protection mechanisms. For a more thorough approach, which involves a deconstruction of various relevant elements and their integration into a multiscale model, see the recent paper of Taylor et al. [10]. Another promising approach, aimed towards the rational design of new corrosion inhibitors, utilizes machine-learning methods (e.g., see the recent review of Winkler [11]), with which it is possible to generate reasonably robust and predictive quantitative models, but the success comes at the expense of much lower mechanistic insight compared to computationally intensive physics-based methods [12,13].
This paper is Part I of a two-part article series (for Part II, see ref [14]) that represents the continuation of the research of our group on the characterization of the adsorption of imidazole, triazole, and tetrazole-used as archetypal models of azole corrosion inhibitors-on copper surfaces by means of DFT calculations as to provide an atomic-scale insight into the chemistry of azole-copper bonding. Molecular structures of imidazole, triazole, and tetrazole are shown in Scheme 1. In previous publications, the adsorption of these molecules was characterized on metallic Cu(111) [15,16] as well as on Cu 2 O(111) [17]. It was shown that neutral molecules bind weakly to Cu(111) via unsaturated N heteroatoms through σ-type bonding and the magnitude of adsorption energy decreases from imidazole to tetrazole. In contrast to neutral molecules, deprotonated molecules bind strongly to Cu(111). This is particularly true for triazole and tetrazole, which form two strong N-Cu bonds, whereas imidazole cannot due to steric reasons associated with its molecular geometry. This observation suggested that, for imidazole, the neutral form and, for triazole and tetrazole, their deprotonated forms are the active species for inhibiting corrosion [16]. Experimental measurements partly supported this inference [18]. However, oxide-free copper surfaces are more relevant at acidic pH, but, under other conditions, copper surfaces are often oxidized [19]. To this end, we have chosen Cu 2 O as a model of the oxidized copper surface and, in our previous publication [17], we characterized the non-dissociative adsorption of these molecules on Cu 2 O(111) and Cu 2 O(111)-w/o-Cu CUS surfaces. In the current two-part article series, we extend that study in two ways: first, we consider additional Cu 2 O surfaces and, secondly, we also consider the feasibility of dissociative adsorption that results from the cleavage of the N1-H bond upon adsorption; this mode of adsorption was not considered in our previous publication.
The following Cu 2 O surfaces are considered within the current two-part article series:  Figures 1 and 2. Among them, only the surface model (i) is stoichiometric, while the rest are non-stoichiometric. The choice of the currently utilized surface models was motivated by considering stability issues of pristine and molecularly covered surfaces. The stability of various Cu 2 O surfaces in oxygen atmosphere was characterized in detail by Soon et al. [20,21] by means of DFT-GGA calculations. The Cu 2 O(111)-w/o-Cu CUS and Cu 2 O(110):CuO surfaces were found to be the stablest; the former under oxygen-lean and the latter under oxygen-rich conditions. The scanning tunneling microscopy (STM) study of Önsten et al. [22] seems to have identified the Cu 2 O(111)-w/o-Cu CUS structure (according to their nomenclature, this surface was labeled as model B of the (1 × 1) terminated surface). In addition, they also observed the ( √ 3 × √ 3)R30 • reconstructed Cu 2 O(111) surface, which contains one-third of a monolayer of ordered oxygen vacancies, and proposed two surface models: model A, which contains only oxygen-vacancies, and model B, which contains also copper vacancies (this model is currently designated as Cu 2 O(111)-recon-( √ 3 × √ 3)R30 • ). On the basis of combined experimental-computational study of methanol dehydrogenation on ( √ 3 × √ 3)R30 • reconstructed Cu 2 O(111), Besharat et al. [23] recently argued in favor of model B, but then the same group of authors also showed that adsorption data of SO 2 strongly suggest model A with the corollary that the surface likely consists of mixture of A and B phases [24]. However, it should be noted that adsorbates can alter the stability of surfaces-a clear example is shown in Part II [14]-hence it is questionable whether one can reliably deduce the structure (or stability) of bare surfaces on the basis of adsorption data.
Recently, Nilius et al. [25] drew attention to the well-known drawback of GGA, which underestimates band gaps and fails to correctly reproduce absolute positions of the band-edges. They showed that GGA too strongly favors non-stoichiometric Cu 2 O(111)-w/o-Cu CUS against stoichiometric Cu 2 O(111), i.e., GGA predicted cost of non-stoichiometry is too small. Using the HSE hybrid functional, they showed, in agreement with their STM experiments, that stoichiometric Cu 2 O(111) becomes more stable than Cu 2 O(111)-w/o-Cu CUS under oxygen-lean conditions. Although we currently use the GGA approach, which is susceptible to the aforementioned problem (it should be noted that even GGA+U does not perform well for Cu 2 O [26], see Appendix A), we nevertheless find that the stoichiometric Cu 2 O(111) surface is stabilized by adsorbed species-i.e., energetic deficiency of coordinatively unsaturated Cu ions is compensated by their stronger bonding to adsorbates.
In the current Part I, we focus on the molecule-surface bonding as well as on the feasibility of dissociative adsorption at various sites on Cu 2 O surfaces, including the O vacancy sites, and compare the adsorption bonding of azole molecules as corrosion inhibitors to that of chloride as a prototype corrosion activator. In Part II [14], we address intermolecular lateral interactions and thermodynamic stability of identified adsorption structures by means of two-dimensional phase diagrams. Due to obvious modeling reasons, the adsorption is considered at a solid/vacuum interface, although, in the context of corrosion inhibition, it would be more appropriate to consider adsorption at a solid/water interface.

Computational Details
DFT calculations were performed with the generalized gradient approximation (GGA) of Perdew-Burke-Ernzerhof (PBE) [28]. This choice is motivated on the one hand by the observation that the GGA+U method, which is often used for the description of transition-metal oxides, does not substantially improve the band-gap of Cu 2 O [26,29] and on the other hand by the fact that hybrid functionals (e.g., HSE) are computationally way too expensive for our available resources. In order to have some idea on how the +U correction would affect the results, we provide some comparison between the PBE and PBE+U methods in Appendix A.
The pseudopotential method with ultrasoft pseudopotentials (US-PP) was used [30,31]. The Kohn-Sham orbitals were expanded in a plane-wave basis set up to a kinetic energy cutoff of 30 Ry (240 Ry for the charge density cutoff). Brillouin zone (BZ) integrations were performed employing the special-point technique [32] using Marzari-Vanderbilt cold smearing [33] of 0.01 Ry. All calculations were done using the PWscf code from the QUANTUM ESPRESSO distribution [34,35], whereas visualization and molecular graphics were produced by the XCRYSDEN graphical package [36].

Cu 2 O as a Model of Oxidized Copper Surface
In compliance with our previous publication [17], we used Cu 2 O slabs without a metal support underneath as a model of oxidized copper surfaces. Such models are appropriate for cases where the oxide layer on top of metal is not ultrathin (note that the reactivity of few Å thick oxide films supported on metals can be very different from the reactivity of bulk oxides [37]). The average thickness of Cu 2 O oxide layer formed on Cu immersed in 3 wt. % NaCl solution was estimated experimentally to be about 2.2 ± 0.3 nm [38], which seems sufficient to make the current model adequate.   (Figure 1a,e). Its surface trilayer contains two distinct copper ions: a coordinatively saturated (CSA) and coordinatively unsaturated (CUS), labeled as Cu CSA and Cu CUS , respectively (the Cu CUS ions are colored more yellowish). It also contains two distinct oxygen ions: O up (colored brighter) and O dn (darker), where the superscripts indicate that they are located above (up) and below (dn) the surface Cu layer, respectively (O up is CUS and O dn is CSA). It should be noted that the high-symmetry Cu 2 O(111) is not stable, and it relaxes such that the Cu CUS ions displace laterally toward two adjacent Cu CSA ions if the symmetry is broken (cf. Figure 1a,e). The resulting relaxed surface is labeled as Cu 2 O(111)r, where "r" stands for "relaxed". Figure 1b shows the Cu 2 O(111)-w/o-Cu CUS structure, which is the Cu 2 O(111) that lacks the Cu CUS ions (notation "-w/o-Cu CUS " stands for "Cu 2 O(111) without Cu CUS ions"); Cu CUS -vacancies are indicated by white dashed circles. The O ions bellow Cu CUS -vacancies are named as O sub , where "sub" stands for "subsurface".

Description of Considered Surfaces
The last considered Cu 2 O(111)-type structure is the surface observed experimentally by Önsten et al. [22] and referred to by them as model B of the reconstructed surface due to ordered oxygen vacancies, which can be described with the ( Figure 1c); we label this structure as Cu 2 O(111)-recon-( where "recon" stands for "reconstructed". It can be derived from Cu 2 O(111)-w/o-Cu CUS because it lacks the Cu CUS ions, but, in addition, one third of O up ions are missing. The respective O-vacancies are indicated with red crosses in Figure 1c. Each O-vacancy is surrounded by three Cu CUS ions, labeled as Cu Ovac ("Ovac" stands for oxygen vacancy) to distinguish them from Cu CUS ions of pristine Cu 2 O(111); note that the two CUS type Cu ions have structurally different environments. The Cu Ovac ions therefore always appear in triplets and we use the term "Cu Ovac site" to designate the site composed of these three ions. We do not consider Önsten's model A of the ( √ 3 × √ 3)R30 • reconstructed surface-the difference between the two models is that model B lacks Cu CUS ions-because no new types of Cu and O ions are introduced in model A. It should be noted that non-stoichiometric Cu 2 O surfaces containing Cu-vacancies possess a sizable magnetic moment. Our calculated value (about 1.3 µ B per Cu-vacancy) is in fair agreement with that reported by Li et al. [40]. Despite the sizable magnetic moment, the energy difference between magnetic (spin-polarized) and non-magnetic (spin-unpolarized) total energies is rather small, being about 30 meV per Cu-vacancy. When, in addition to Cu vacancies, O vacancies are also present (e.g., Cu 2 O(111)-recon-( √ 3 × √ 3)R30 • model), then both magnetic moments as well as the difference between magnetic (spin-polarized) and non-magnetic (spin-unpolarized) total energies are reduced. Hence, the results presented herein were obtained with non-magnetic surface calculations.

An Issue with Cu CUS Sites: A Cu 2 O(111) w/o +1Cu CUS Model
As mentioned above, the high-symmetry stoichiometric Cu 2 O(111) is not stable: if the symmetry is broken, then each Cu CUS ion relaxes laterally toward two adjacent Cu CSA ions in the direction pointing to nearby O up ions. The resulting Cu 2 O(111)r is by about 4 meV/Å 2 more stable than its high-symmetry analogue [27]. In Figure 1e, all the Cu CUS ions are depicted to relax in the same direction. However, there are three symmetry equivalent directions along which a given Cu CUS ion may relax; these are 120 • apart from each other and point towards the nearby O up ion. Each Cu CUS ion may therefore relax in one of these three directions, leading to different possible Cu CUS relaxation patterns, which, according to our calculations, display slightly different relaxation energies. This can result in some spurious, though small, contribution to adsorption energy because molecular adsorption may further lower the symmetry and result in different relaxation patterns of the nearby Cu CUS ions. To minimize this potential artifact, the surface model used for the adsorption onto Cu CUS site is modified so that all Cu CUS ions, except those that bond to adsorbed molecules, are removed. The resulting model can be seen as a Cu 2 O(111)-w/o-Cu CUS containing Cu CUS defects, where the number of Cu CUS ions equals the number of adsorbed molecules (Figure 2a). This model is designated as Cu 2 O(111) w/o +1Cu CUS , where the subscript "w/o" is a shorthand for "w/o-Cu CUS " and the suffix "+1Cu CUS " conveys that there is one Cu CUS ion per adsorbed molecule. Note that, at a molecular coverage of one molecule per (1 × 1) unit-cell, the Cu 2 O(111) w/o +1Cu CUS model becomes synonymous with the stoichiometric Cu 2 O(111).
The choice of the Cu 2 O(111) w/o +1Cu CUS model can be further justified by the fact that the relaxed stoichiometric Cu 2 O(111)r, though being a local minimum structure, is not thermodynamically stable [27]. However, in our previous study [17], we showed that molecular adsorption is able to compensate the thermodynamic deficiency of Cu CUS ions. Hence, it seems appropriate to maintain only those Cu CUS ions that bond with adsorbed molecules.
Cu CUS ions display one further interesting characteristic: on the bare surface, they relax laterally as discussed above (cf. Figure 1e). However, when a Cu CUS ion bonds with an adsorbed molecule, then it shifts back to the high-symmetry position ( Figure 2b). Apparently, the adsorbed molecule sufficiently diminishes the unsaturated character of Cu CUS , thus withdrawing its need to bond with nearby Cu CSA ions. To better reflect the molecule-surface bond strength, the binding and adsorption energies at Cu CUS sites are calculated with respect to the high-symmetry position of Cu CUS ion on the pristine surface. The energy difference between high-and low-symmetry position of Cu CUS ion is 0.13 eV for Cu 2 O(111)-(1 × 1), 0.04 eV for Cu 2 O(111) w/o +1Cu CUS -(2 × 2), and vanishes for

Molecular Labels
The labels MolH, Mol − , and Mol • designate azole molecules in neutral, deprotonated, and radical form, respectively, whereas Mol is used as a generic label for a dissociated molecule-a molecule with the N1-H bond broken and stripped of the pertinent H-when the charge of the species is not of concern. ImiH, TriH, and TetH are shorthand labels for intact imidazole, triazole, and tetrazole, respectively. For other speciation forms, analogous designations as for MolH are used (cf. Scheme 1), e.g., for imidazole, the respective labels are Imi − , Imi • , and Imi.  [41] was applied to cancel an artificial electric field that develops along the direction normal to the slab due to periodic boundary conditions imposed on the electrostatic potential. Molecular adsorption was modeled with the (2 × 2) supercell for Cu 2 O(111) based models, (
We consider molecular (or non-dissociative) as well as dissociative adsorption of azole molecules. In the latter, the N1-H bond (cf. Scheme 1) is broken upon adsorption. Non-dissociative adsorption of azole molecules is described as: whereas dissociative adsorption can be written as: where standalone * designates a free adsorption site, while MolH * , Mol * , and H * denote adsorbed species. The respective molecular and dissociative adsorption energies (per adsorbed molecule) are calculated as: and where E MolH , E slab , E MolH/slab , and E Mol+H/slab are the total energies of isolated intact MolH molecule, clean slab, MolH/slab, and coadsorbed Mol+H/slab system, respectively. To estimate how strong the dissociated molecule, H, and Cl bind to the surface, we utilize the binding energy (E b ) or the bond-strength (D); the E b is calculated as: where A stands for adsorbate (Mol, H, or Cl) and the meaning of the energy terms is analogous to those defined above. By convention, the bond-strength D is positive and hence opposite to binding energy, D = −E b . Note that E b is calculated with respect to isolated radicals (Mol • , H • , or Cl • ) and not ions (Mol − , H + , or Cl − ); isolated radical species were calculated with spin-polarized calculations.
Note that for non-dissociative adsorption the binding energy is equivalent to adsorption energy, i.e., The relative stability of molecular vs. dissociative adsorption can be evaluated by considering the dissociation reaction occurring at the surface: The respective dissociation energy is calculated as: and dissociative adsorption is favored over the non-dissociative adsorption when ∆E < 0.

Charge Density Difference
The formation of a chemical bond between the adsorbate and the surface can be characterized in terms of the charge density difference, ∆ρ(r), defined as: where the subscripts have the same meaning as in Equation (5). For ∆ρ(r) calculations, the geometries of the "slab" and standalone "A" structures are kept the same as in the "A/slab" system.

Bader Charge Analysis
Bader charges [42] were calculated by first generating charge densities with single point selfconsistent-field calculations of US-PP optimized structures using the projector-augmented-wave (PAW) potentials [43] with 40 Ry and 1000 Ry kinetic energy cutoffs for wave-functions and charge densities, respectively, and then computing the Bader charges using the bader program [44,45].

Dissociation Activation Energies
Dissociation activation energies were calculated using the climbing-image nudged-elastic-band (CI-NEB) method [46,47] that models an elementary reaction step as the minimum energy path (MEP) connecting the initial state (IS) with the final state (FS). The configuration with the maximum energy along the MEP is the transition state (TS) and the activation energy is calculated as the difference between the TS and IS energies, E TS − E IS . For the precise location of the TS, the threshold for the magnitude of the atomic forces was set below 40 meV/Å.

Adsorption of Standalone Molecules at Lower Coverage
The focus of this section is on the molecule-surface interactions (lateral intermolecular interactions will be considered in Part II [14]), hence we will analyze the adsorption at lower coverage. By lower coverage, we do not mean an extra low coverage, but instead some intermediate coverage, which is low enough that it allows for focusing predominantly on the molecule-surface interactions, even though at such coverage the long-ranged dipole-dipole interactions [15,17,48] may not yet be insignificant. To this end, calculations were performed by adsorbing a single molecule in the following surface supercells: (2 × 2) for the two Cu 2 O(111) based structures (Figure 3a   Intact molecules (MolH) typically adsorb with the molecular plane perpendicular to the surface (or nearly so) and form one N-Cu bond and one X-H· · · O hydrogen bond (X = C2 for ImiH or N1 for TriH and TetH), except when bonding to Cu Ovac sites (Figure 6), where ImiH binds with the N3 atom to two adjacent Cu Ovac ions and forms the C2-H· · · O sub hydrogen bond. TriH and TetH would also adsorb analogously via the N2 atom to two adjacent Cu Ovac ions and form the N1-H· · · O sub hydrogen bond; however, in this geometry, the N1-H bond breaks barrierlessly for TriH * and TetH * (this issue is considered in more detail in Section 3.3.3). Hence, alternative adsorption forms, which are less susceptible to N1-H bond cleavage, are depicted in Figure 6; in these forms, the TriH * and TetH * bind in a tilted geometry without hydrogen bonding and form two N-Cu Ovac bonds, such that each of the two N atoms (N2 and N3) binds to a different Cu Ovac ion.        It is evident from Figure 8a that all three molecules display similar non-dissociative adsorption energies; the largest binding energy difference between the two intact molecules adsorbed at a specific site is displayed by the ImiH and TetH bonded to Cu Ovac , 0.25 eV. Figure 8 further reveals that, in general, the molecules bind stronger to CUS than to CSA sites, although the specific geometry of the adsorption site matters, i.e., the two CUS and the two CSA sites display different reactivity and the trend of the non-dissociative adsorption bond strength is: Cu CUS Cu Ovac > Cu surf > Cu CSA . In particular, the bonding of intact molecules to Cu CUS sites is rather strong with the bond strength of about 1.5 eV, whereas at other sites it is considerably weaker, between 0.7 and 1 eV at Cu Ovac , about 0.7 eV at Cu surf , and about 0.5 eV at Cu CSA . Strong molecular adsorption at Cu CUS site is evident also from the N-Cu bond lengths, which are about 1.9 Å at Cu CUS and more than 2.0 Å at other sites.

Cu 2 O(111)-w/o-Cu
In contrast to non-dissociative adsorption modes, where all the molecules display similar binding energies, considerable differences between the molecules appear for dissociative adsorption modes. In particular, dissociated imidazole (Imi * ) binds considerably weaker than dissociated triazole (Tri * ) and tetrazole (Tet * ). The reason is due to the molecular geometry: triazole and tetrazole have adjacent N atoms, but imidazole does not. Hence, Tri * and Tet * bond with at least two N atoms to Cu sites, whereas Imi * usually binds to the surface only with a single N atom because the other N atom is located on the other side of the molecule and is not available for bonding with the surface when imidazole adsorbs with its molecular plane perpendicularly to the surface (or nearly so). Similar differences between imidazole and triazole/tetrazole were also observed on metallic Cu(111) [16]. Note that in some cases Imi * can bond with both N atoms to the surface, provided its molecular plane is sufficiently tilted, but the resulting two N-Cu bonds are frustrated and such bonding is only by about 0.2 eV stronger than the single-atom perpendicular boning; e.g., this happens on Cu 2 O(111)-w/o-Cu CUS ( Figure 5) and on Cu(111) [18]. The single N atom bonding of imidazole is also the reason that the difference in adsorption bond-strength between dissociated Imi * and intact ImiH * is considerably smaller compared to those displayed by triazole and tetrazole, i.e., for imidazole the Mol * vs. MolH * the difference ranges from 0.1 eV at Cu surf to 1.1 eV at Cu Ovac , whereas, for triazole and tetrazole, it ranges from about 1.6 eV at Cu CUS to about 2.9 eV at Cu Ovac . The difference between the MolH * and Mol * binding is the largest at the Cu Ovac site because the Mol * species bind the strongest to it, but the MolH * species bind the strongest to Cu CUS instead. Namely, for intact MolH species, a single N-Cu bond is sufficient for effective adsorption (note that Cu CUS ion is standalone), whereas, for strong adsorption of dissociated Mol * species, at least two N-Cu bonds are required. The Cu Ovac site is therefore preferred for Mol * species because it consists of three adjacent Cu Ovac ions (see Figure 1) and can form two strong N-Cu Ovac bonds (Tri * and Tet * form one single and one bifurcated N-Cu Ovac bond, Figure 6), but, at the Cu CUS site, the Mol * species can form only one N-Cu CUS bond; Tri * and Tet * thus form the second bond with the adjacent Cu CSA ion (Figure 4), which is considerably less reactive. The trend of the adsorption bond-strength of Mol * at different surface sites is therefore Cu Ovac > Cu CUS Cu CSA > Cu surf (Figure 8). For Imi * , the respective binding energies are about −2.1, −1.9, −0.9, and −0.7 eV at Cu Ovac , Cu CUS , Cu CSA , and Cu surf , respectively, whereas, for Tri * and Tet * , these values are about −3.7, −3.1, −2.2, and −2.0 eV, respectively.  Figure 9 with the bottom part of Figure 6).

Electronic Structure Analysis
To shed some light onto the nature of the molecule-surface bonding, we utilize the charge density difference, ∆ρ(r) of Equation (9), which is presented in Figure 10, where red (blue) color represents the electron charge excess (deficient) regions; top-row shows the MolH-surface and bottom-row the Mol-surface bonding. For the sake of brevity, the ∆ρ(r) is plotted only for the triazole molecule, which was chosen because its electronegativity and chemical-hardness are in between that of imidazole and tetrazole [15] (triazole can thus be seen to represent the average behavior of the three).
As for the TriH-surface bonding (top-row), the ∆ρ(r) plots clearly reveal the N-Cu bonds (note the red electron charge accumulation lobe in the midst of the N-Cu bonds) as well as the N1-H· · · O hydrogen bonds (note the red charge accumulation region above the O ion in the direction towards the H atom of the triazole). The intensity of the red charge accumulation lobe in the midst of the N-Cu bonds follows the reactivity trend of Cu ions towards the bonding to MolH, i.e., Cu CUS > Cu Ovac > Cu surf > Cu CSA . Indeed, the N-Cu CSA bond clearly displays the faintest red lobe. Bader population analysis reveals that the charge of TriH * is slightly positive on all considered Cu sites, ranging from +0.03 on Cu 2 O(111)-recon-( √ 3 × √ 3)R30 • site to +0.08 on Cu 2 O(110):CuO. As for dissociated adsorption modes, the ∆ρ(r) plots reveal that Tri-surface bonding is stronger than the TriH-surface bonding because the charge accumulation and deficit regions are more intense for the former, in particular, the red electron charge accumulation lobes in the midst of the N-Cu bonds are more intense for Tri * compared to TriH * . According to Bader population analysis, the Tri * is negatively charged, but it is not fully anionic. Instead, it is about midway between the Tri • radical and Tri − anion. In particular, molecular charges are −0.58, −0.56, −0.65, and −0.47 for Tri * bonded to Cu CUS , Cu CSA , Cu Ovac , and Cu surf , respectively. These Bader charges moderately correlate with Tri * binding energies at these sites (the correlation coefficient is 0.91), i.e., the larger the negative charge of Tri * is, the stronger is its bonding. The charges obtained on the first three sites, which belong to Cu 2 O(111) type surfaces, are similar to the Tri * charge of −0.61 obtained on Cu(111) [16]. .006 e/a 3 0 and in (g) also the ±0.006 e/a 3 0 isosurfaces are shown. The blue (red) color represents the electron deficient (excess) regions, i.e., electron charge moved from blue to red regions. For Tri * on Cu 2 O(111)-w/o-Cu CUS , the second stablest identified adsorption structure is shown instead of the stablest because the latter is highly tilted and forms three out-of-plane N-Cu bonds (hence, it is not possible to plot the three bonds simultaneously on a single contour plotting plane).
We further calculated Bader charges for all the three molecules adsorbed on the Cu CUS site (Table 1). For intact MolH * , Bader charges are +0.13, +0.05, and +0.02 for imidazole, triazole, and tetrazole, respectively, whereas, for dissociated Mol * , Bader charges are −0.41, −0.58, and −0.64 for Imi * , Tri * , and Tet * , respectively. Bader charges thus follow the imidazole > triazole > tetrazole trend (i.e., tetrazole has the most negative charge) for both MolH * and Mol * adsorption modes, which is in accordance with their Mulliken electronegativity, i.e., imidazole is the least and tetrazole the most electronegative [16]. Respective correlations between molecular Bader charges of adsorbed molecules and their electronegativities are shown in Figure 11a (Mulliken electronegativities are taken from Ref. [16] and were calculated from vertical ionization potentials (I) and electron affinities (A) that were obtained from the X → X + + e − and X − → X + e − reactions, respectively, where X ≡ Mol or MolH). It is worth noting that these molecular Bader charges show no relation to the adsorption bonding trend for the three molecules, which may seem surprising although such behavior is not unknown. For example, Stenlid et al. [49] observed recently in their study of several probe molecules on TiO 2 nanoparticles that larger charge transfer does not necessarily lead to stronger interaction energies. In contrast to molecular charges, Bader charges of bare Cu ions before molecular adsorption (+0.72, +0.68, +0.34, and +0.29 for Cu surf , Cu CSA , Cu CUS , and Cu Ovac , respectively) correlate remarkably well with adsorption binding energies of dissociated molecules, i.e., the smaller the charge of the Cu ion is, the stronger is its bonding with the Mol * (Figure 11b). Dissociated molecules are negatively charged in the adsorbed state and it seems that the less the Cu ion is positively charged, the easier is the charge flow to the molecule and consequently the Cu ion is more susceptible to bonding with Mol * .
In contrast, such correlations are much weaker for neutral molecules-the correlation coefficients are only 0.80, 0.67, and 0.62 for ImiH * , TriH * , and TetH * , respectively-and the reason for this can be associated with the following two observations: (i) the adsorption induced charge transfer for MolH * is very small and in the opposite direction with respect to that of Mol * and (ii) the MolH-surface bonding also involves an H-bond, which is due to the availability of the nearby O ion, and it can be seen from the top row of Figure 10 that some sites have a more appropriate geometry than others for the formation of an H-bond. Although the adsorption energies of MolH * and Bader charges of bare Cu ions do not show the same trend-the adsorption energies display the Cu CUS < Cu Ovac < Cu surf < Cu CSA and Cu Bader charges the Cu Ovac < Cu CUS < Cu CSA < Cu surf trend-it is possible to make a classification into the following two groups: (i) unsaturated Cu CUS and Cu Ovac ions are less positively charged and bind MolH stronger, whereas (ii) saturated Cu surf and Cu CSA ions are more positively charged and bind MolH weaker.

Dissociation of Adsorbed Azole Molecules, MolH * + * → Mol * + H *
Despite the fact that dissociated Mol * species bind considerably stronger to the surface than intact MolH * species (Figure 8), the stronger bonding is not sufficient to compensate for the cost of the N1-H bond cleavage in the majority of currently considered cases. This is evident from Figure 12, which plots the dissociation energy for the MolH * + * → Mol * + H * reaction. The respective dissociation reaction is exothermic only if E diss ads < E ads , cf. Equation (8), and this condition is met only for triazole and tetrazole bonded to O vacancies, that is, Cu Ovac sites on Cu 2 O(111)-recon-( √ 3 × √ 3)R30 • and Cu Ovac (110) sites on Cu 2 O(110):CuO-w-Ovac (the latter sites are not considered in Figure 12). The results for triazole are in good agreement with those published previously for benzotriazole [27]: the ∆E values for both molecules are about 0.1, 0.4, and 0.5 eV at Cu CSA , Cu CUS , and Cu surf sites, respectively. Figure 12 further reveals that dissociation is by far the most endothermic for imidazole, which is due to the small difference in adsorption bond-strength between intact ImiH * and dissociated Imi * discussed above; this small difference stems from single N atom bonding of Imi * .

Bonding of H to Various Sites on Cu 2 O Surfaces
To complete the description of the dissociation reaction, MolH * + * → Mol * + H * , we also need to describe the bonding of H to various sites on Cu 2 O surfaces. It seems intuitive to assume that H * bonds to O surface sites. Calculations indeed confirm this anticipation; however, they also reveal that H can bond strongly to unsaturated Cu ions (Table 2), in agreement with previous studies [23,50]. Notably, H binds the strongest to hollow site consisting of three adjacent Cu Ovac or Cu Ovac  Table 2 summarizes the calculated adsorption data for H adsorbed at various sites on Cu 2 O surfaces, whereas the corresponding structures are depicted in Figure 13. Our calculated H-O and H-Cu bond lengths are in good agreement with those reported in the literature [23,50]. Bader analysis reveals that although H bonded to O sites is positively charged, it is not fully a proton-its Bader charge is about +0.6. In agreement with Yu et al. [50], we also find that the Bader charge of H adsorbed to unsaturated Cu sites is significantly negative, being about −0.2 at Cu CUS and −0.3 at Cu Ovac and Cu Ovac (110) , while the charge of H adsorbed at saturated Cu sites is close to zero.

Co-Adsorption of Mol and H
Comparison of Mol * data presented in Figure 8 and H * data presented in Table 2 clearly reveals that, after the dissociation, the preferred co-adsorption structures consist of Mol * bonded to Cu sites and H * bonded to O sites; the latter thus forms an OH group at the surface. This is true even for Cu Ovac and Cu Ovac (110) sites, to which H * binds stronger than to O sites. However, the H * preference for Cu Ovac sites (0.3 eV) is considerably smaller compared to that of Mol * (about 1.4 eV; this number is the difference between the E b of Mol * on Cu CSA and Cu Ovac sites). The dissociation of adsorbed azole molecules can therefore be written such that the formation of an OH * group is indicated: where O * indicates a given lattice O ion at the surface of Cu 2 O. Figure 14 schematically illustrates the mechanism of dissociative adsorption of triazole and tetrazole on oxygen vacancy Cu Ovac and Cu Ovac (110) sites along with the involved energy differences. These two sites are specifically considered because the dissociation is favorable only thereon (cf. Figure 12). According to our calculations, the dissociation of triazole and tetrazole on these sites is very facile. Namely, for a properly oriented molecule, the dissociation of the N1-H bond proceeds either without a barrier or with a marginal energy barrier ( Figure 15). In particular, the N1-H bond cleavage is barrierless for triazole on Cu Ovac site and for tetrazole on both types of oxygen vacancy site, whereas, for triazole on Cu Ovac (110) , the dissociation barrier is only 11 meV. The N1-H bond cleavage is therefore much easier at oxygen vacancy sites on Cu 2 O compared to pristine metallic Cu surfaces, where the energy barrier was calculated to be about 1 eV for benzotriazole [51]. It is worth noting that oxygen vacancy sites were found to be reactive also for the cleavage of the N-H bond of benzotriazole [52] and the O-H bond of methanol [23]. However, at few other sites on Cu 2 O surfaces, where dissociation of azole molecules is endothermic, the opposite occurs, i.e., on Cu CSA and Cu CUS sites a properly oriented Mol * + OH * transforms into MolH * + O * either barrierlessly for triazole or with a vanishing barrier of about 10 meV for tetrazole.    Figure 15. Beware that the two shown TetH * structures are not local minima and were obtained by constraining the N1-H bond-length; also triazole at Cu Ovac is not a local minimum on the PES but instead a wide plateau. These three structures are shown as faded to indicate their instability. The shown TriH @ Cu Ovac and TetH @ Cu Ovac structures display more exothermic adsorption energies than those considered in Figure 6, but note that the current structures are not stable minima on the PES.

Mechanistic Insight into MolH * + O * → Mol * + OH * Dissociation
Finally, it should be noted that Bader analyses (Tables 1 and 2) reveal that it is not fully appropriate to designate the dissociation of adsorbed azole molecules as the MolH ads → Mol − ads + H + ads because the charge of the two product species is significantly different from −1 and +1. This is why we write the dissociation of adsorbed azole molecule either as MolH * + * → Mol * + H * or as MolH * + O * → Mol * + OH * , where the suffix * indicates the adsorbed species without any reference to its charge. TriH @ Cu Ovac TetH @ Cu Ovac Figure 15. Calculated minimum energy paths on the PES for the dissociation of TriH * and TetH * on oxygen vacancy sites, i.e., the elementary step #2 in Figure 14. Beware that due to small energy differences the energy unit of the ordinate axis is meV and not eV. The labels (A-D) comply with the labels written below the MolH * structures in Figure 14.

Comparison between Molecular and Chloride Adsorption
The comparison between azole and chloride adsorption is of interest because corrosion is usually promoted by some reactive corrosive species and chloride can be seen as a prototypical corrosion activator. Moreover, the inhibitive effect of azoles for corrosion of copper in chloride media has often been investigated (e.g., see ref [2]). Figure 16 shows the most stable identified adsorption structures of Cl * at various Cu sites on Cu 2 O surfaces along with the respective binding energies and Bader charges. Similarly as reported above for Mol * , Cl * is also not fully anionic. Instead, its Bader charge is about −0.5 at all considered sites. Cl * binds the strongest to oxygen vacancy Cu Ovac and Cu Ovac (110) sites (E b = −3.8 eV), followed by the Cu CUS sites (E b = −3.0 eV), whereas the bonding to saturated Cu CSA and Cu surf ions is considerably weaker (E b = −2.0 and −1.8 eV, respectively). These Cl-surface bond strengths are considerably stronger than those displayed by intact azole molecules (cf. Figure 8a). Only triazole and tetrazole molecules in dissociated Mol * form interact with the surface strongly enough to rival the Cl-surface bonds. A similar trend was also observed on metallic Cu surfaces [6,53]. To facilitate the comparison between the Cl * and Mol * adsorption bonding, Table 3 tabulates the respective binding energies. Note that the bonding of Imi * is considerably weaker compared to that of Cl * , Tri * , and Tet * and the reason was already explained above, i.e., for strong bonding, the Mol * has to form at least two strong N-Cu bonds, which is not possible for imidazole due to its incompatible geometry. As for the comparison of E b between Cl * , Tri * , and Tet * , the molecules bond somewhat stronger than Cl * to Cu CSA , Cu surf , and Cu CUS sites, whereas the opposite is true on Cu Ovac and Cu Ovac (110) sites. The adsorption bonding of Tri * , Tet * , and Cl * is therefore comparable in strength and this may be of relevance for the competitive adsorption scenario as a plausible mechanism of corrosion inhibition, i.e., Tri * and Tet * may hinder the adsorption of Cl * ; note that chloride is known to induce faster thinning and eventual breakdown of the passive film, followed by pit nucleation [54,55]. However, to go beyond this crude qualitative statement is not appropriate because current calculations were performed in a vacuum, while corrosion typically occurs at the solid/water interface.

Adsorption of Molecules from Vacuum and Aqueous-Phase: Differences
The above comparison between molecular and chloride adsorption reveals that only Tri * and Tet * adsorb strong enough to rival the Cl-surface bonds. However, there is an important point to keep in mind when considering the results presented above in the context of corrosion inhibition. Namely, due to obvious modeling reasons, the calculations were performed at a solid/vacuum interface, whereas, for corrosion, the solid/water interface is relevant. The adsorption at the latter interface is competitive (or substitutional) because the surface is always covered with solvent molecules and also with other species, such as hydroxyls. Thus, a given molecule will adsorb only if its adsorption is competitive enough to substitute other species from the surface. In contrast, at the solid/vacuum interface, the surface is clean and the molecule adsorbs readily unless its interaction with the surface is repulsive. Furthermore, the solvent considerably affects the energetics of adsorption because during the adsorption the molecule must get rid, at least partially, of its solvation shell and also displace solvent molecules from the surface. To get the first idea about the effect of solvent on adsorption energetics, one can compare the molecular solvation free energy with the in vacuo adsorption energy: if the former is significantly stronger than the latter, then the molecule is unlikely to adsorb, because it interacts more strongly with the solvent than the surface. The solvation free energies, as calculated by the COSMO (conductor-like screening model) implicit solvent model, are about −3 eV for Mol − and Cl − and about −0.5 eV for MolH [16,53]. Comparing these values with the corresponding binding energies implies that Mol − and Cl − are unlikely to adsorb at coordinatively saturated Cu sites on Cu 2 O surfaces not only because they instead prefer to bond to unsaturated Cu sites, but moreover because these species interact more strongly with the solvent than with the CSA sites (i.e., solvation free energies are about −3 eV, while CSA binding energies are about −2 eV).
Current results show that dissociated Mol * species bond considerably stronger to Cu 2 O surfaces than the intact MolH * molecules, which is due to the more reactive nature of Mol species that stems from its dangling bond at N1. This argument is not hindered by the fact that the binding energies were calculated with respect to the isolated Mol • radical (see Appendix B). However, despite the much stronger adsorption bonding of Mol * and Cl * species, we can reasonably anticipate that the net adsorption energy difference between MolH and Mol − is considerably diminished in aqueous-phase [16,52] because charged Mol − and Cl − species interact strongly not only with the substrate, but also with the solvent. This implies that adsorption from the aqueous-phase involves the change from one stable environment (solvent) to another one (surface), hence the adsorption energy is to a significant extent given by the net difference between the two strong interactions. We can therefore reasonably infer that the net aqueous-phase adsorption energies of Mol − and Cl − are considerably smaller in magnitude than the respective in vacuo energies. , respectively. We find that, in general, unsaturated Cu sites bond adsorbates much stronger than saturated sites. Among them, the O vacancy sites, which consists of a triplet of unsaturated Cu ions, are found particularly reactive towards the dissociation of triazole and tetrazole (N1-H bond cleavage); these sites display similar characteristics on both (111) and (110) surfaces, which we attribute to their similar local geometry. The dissociation of triazole and tetrazole at O vacancy site is assisted by the adjacent surface O ion because the N1-H bond cleavage proceeds by H transfer from the molecule to the nearby O ion, thus forming a surface OH group: the bond cleavage proceeds barrierlessly (or nearly so) thus being remarkably easier than on pristine metallic Cu surfaces, where the energy barrier is on the order of 1 eV for benzotriazole. Despite this large difference in a dissociation activation barrier, Bader analysis reveals that dissociated triazole and tetrazole display a similar charge on Cu 2 O and on pristine Cu(111) surfaces, about −0.6. Indeed, Bader charges of adsorbed molecules show no relation to the strength of the molecule-surface bond. Instead, the Bader charges of bare Cu ions before molecular adsorption correlate remarkably well with the adsorption binding energies of the dissociated molecules, i.e., the smaller the charge of the Cu ion, the stronger is the molecule-surface bond.

Conclusions
While all three azole molecules display similar non-dissociative adsorption energies, dissociated triazole and tetrazole adsorb considerably stronger than dissociated imidazole, which is also why the dissociation is favorable only for triazole and tetrazole. Indeed, dissociated triazole and tetrazole adsorb strong enough to rival the Cl-surface bonding. While it is tempting to associate this dissociation tendency of the three molecules to their pKa constants (note that imidazole is considerably more basic than triazole and tetrazole and consequently less susceptible to deprotonation; their pKa constants at 25 • C for MolH (aq) Mol − (aq) + H + (aq) are 14.5 [56], 9.4 [57], and 4.7 [58], respectively), we believe that this correlation is coincident in the current case. This tendency is instead related to molecular geometry because triazole and tetrazole can form two strong N-Cu bonds, but imidazole cannot due to its incompatible molecular geometry. This observation may be of relevance to corrosion inhibition studies, when considering the speciation of a given inhibitor-i.e., which speciation form is in the majority at a given pH-because the active molecular species at the surface may differ from the molecular form that is in the majority in the solution not only because the pH in the bulk solution and the local pH near the surface can differ significantly, but also because the chemical characteristics of the solvent and the surface are different and, moreover, because the solvent is a 3D system while the surface is a 2D system. The latter implies that geometric steric requirements of the two can differ appreciably and, precisely in this aspect, the imidazole differs from triazole and tetrazole.
Author Contributions: Anton Kokalj and Dunja Gustinčič conceived and designed the computational study. Dunja Gustinčič performed a large majority of the calculations, while Anton Kokalj performed NEB, GGA+U, and magnetic surface calculations. Data were analyzed and the paper was written jointly by both authors.

Acknowledgments:
The authors thank Matic Lozinšek and Matic Poberžnik for valuable discussions and for their careful reading of the manuscript. The authors acknowledge the financial support from the Slovenian Research Agency (Grant No. P2-0393). This work is a part of the M-Era.Net project "COIN DESC: Corrosion inhibition and dealloying descriptors", co-financed by the Slovenian Ministry of Education, Science, and Sport (Grant No. C3330-17-500074).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Effect of Hubbard U Parameter on Electronic and Adsorption Properties
Given that the results presented in this paper are obtained with plain GGA, which is known to underestimate band gaps and fails to correctly reproduce absolute positions of band-edges, we present below the effect of the Hubbard U parameter on the electronic and adsorption properties. To this end, we have used the simplified version of GGA+U method of Cococcioni and de Gironcoli [59]. The GGA+U method is often used, due to its computational efficiency, to correct for the aforementioned problems, although, for copper-oxides, it does not perform well [26]. Figure A1a plots the band gap of Cu 2 O bulk as a function of U eff = U − J parameter (applied to Cu d orbitals), where U and J are parameters describing screened on-site Coulomb and exchange interactions, respectively. Note that the values used in the literature for the U eff parameter of Cu 2 O range from 3 to 8 eV [23,24,26,29,60,61]; Yu and Carter [60] reported the ab initio calculated value of 3.6 eV. It is evident from the figure that the +U correction only marginally increases the band gap, i.e., from the value of 0.43 eV for U eff = 0 eV (plain PBE) to the value of 0.68 eV at U eff = 9 eV, in fair agreement with previous literature reports [26,29], whereas the experimental value is 2.17 eV [62]. Although GGA+U fails to significantly correct the band gap problem, further inspection reveals ( Figure A2) that it downshifts the position of the Cu d-band. As for the structural properties, Figure A1b shows that the Cu 2 O lattice parameter increases with increasing U eff , from 4.34 Å for plain PBE to 4.39 Å for U eff = 9 eV. According to our calculations, the +U correction thus slightly worsens the agreement with experiment because the experimental value is 4.27 Å. This is in variance with previous studies [26,29,60,61], which reported that the lattice parameter decreases with the increasing value of U eff . We therefore made several further tests (also using the PAW potentials), but the lattice parameter always increased with the increasing U eff .  Finally, we address the effect of the +U correction on the adsorption characteristics of triazole at three different sites on Cu 2 O(111) type surfaces. Given that our lattice parameter increases with U eff , which is in variance with literature reports, we also tested the effect of lattice parameter on the adsorption properties. Hence, we made calculations for two sets of lattice parameters, i.e., at the PBE+U optimized values and at the value given by the PBE. Differences between the two sets of adsorption results were insignificant and we report below results only for the GGA+U optimized lattice parameters. Figure A3a plots non-dissociative adsorption energies (E ads ) of triazole at Cu CSA , Cu CUS , and Cu Ovac sites. It can be seen that the +U correction has almost no effect on the adsorption energy at the least reactive Cu CSA site, whereas at more reactive unsaturated Cu CUS and Cu Ovac sites, the adsorption bond strength increases with increasing U eff . In particular, it increases by 0.1 eV and 0.2 eV for Cu CUS and Cu Ovac sites, respectively, as passing from U eff = 0 to 6 eV. In addition to the enhancement of the molecule-surface bond, the +U correction also favors the dissociation of triazole at these sites ( Figure A3b), i.e., the exothermicity of dissociation increases by about 0.2 eV as passing from U eff = 0 to 6 eV (note also that at about U eff = 3 eV dissociation at Cu CSA becomes favorable). Apart from these differences, Figure A3 further suggests that the relative stability of molecular adsorption at various sites is not affected by the +U correction.
The bottom line of this analysis is that the +U correction does not remedy the band-gap problem of Cu 2 O bulk, slightly changes its lattice parameter, moderately enhances the molecule-surface bond strength, relatively favors molecular dissociation, and finally does not alter the reactivity trend of the adsorption sites. Given that one of the principle objectives of this paper is to emphasize the importance of dissociated molecular forms, in this case, using GGA can be considered a more conservative approach because GGA+U relatively favors dissociation with respect to GGA.  Figure A3. Effect of Hubbard U eff parameter (applied to Cu d orbitals) on (a) non-dissociative adsorption energy (E ads ) of triazole and (b) dissociation energy (∆E) of triazole at Cu CSA , Cu CUS , and Cu Ovac sites. Calculations for Cu CSA and Cu CUS sites were performed with one molecule per (2 × 2) supercell, whereas, for the Cu Ovac site, the ( √ 3 × √ 3)R30 • supercell was used, i.e., molecular coverages are the same as those shown in Figure 3.