Density Functional Theory: An Essential Partner in the Integrated Computational Materials Engineering Approach to Corrosion

Computational modeling is increasingly playing a role in advanced approaches to materials design, characterization, lifetime prediction, and risk assessment. Advances in the speed and versatility of modeling algorithms, combined with continued development of computer resources, ensure that the role of computation will continue to grow. The expansion of algorithmic and resource capabilities has already resulted in changes of engineering practices and institutional policies. Take, for example, the use of CALculation of PHAse Diagram (CALPHAD) to generate phase diagrams based on thermodynamic databases catalogued from either experiment or the use of high level ab initio calculations.¹  More generally, the fields of integrated computational materials engineering (ICME) and the Office of Science and Technology Policy’s Materials Genome Initiative point the way toward a widespread integration of computational modeling and data sciences with traditional materials science approaches in the very near future.^2-5  Advanced computational models have also been applied to corrosion science and engineering, ranging from the ab initio modeling of corrosion inhibitors through to the probabilistic risk assessment of pipeline corrosion via Monte Carlo simulation and Bayesian network analysis.^6-7  One tool that has started to be applied more frequently for characterizing corrosion mechanisms is density functional theory (DFT).⁸  DFT is a first principles method capable of evaluating the physical nature of any atomic configuration or transformation and, as such, will continue to play an important role, whether the goal be a well-defined mechanistic study or an integrated multiscale model. Due to the increasing application of this method in the corrosion literature, it seems an appropriate time to provide a brief introduction to the method, a discussion of the kinds of systems that can be approximated in a DFT model, and some case studies in which this method has been usefully applied. This paper also provides a look forward to future opportunities in which DFT can be anticipated to impact the fields of corrosion science and engineering.

DFT is an electronic structure method. Its usefulness to problems in the chemical and materials sciences derives from the fact that the thermodynamic ground states of both molecules and materials consist of a distribution of electrons in equilibrium with a field of positively charged point nuclei. Because the mass of a nuclei is much higher than an electron, the time scale for electronic relaxation is many orders of magnitude faster than the motion of nucleus. Therefore, it can be assumed that the calculation of motion of atomic nuclei and electrons can be separated. This picture is called the Born-Oppenheimer approximation.⁹  When the positions of the nuclei move about, the electronic distribution also shifts in such a way as to minimize the total system energy. Methods that can simulate the electronic structure of molecules or materials, therefore, should be considered among the most rigorous of first principles techniques to predict the thermodynamic and kinetic properties of matter.

Monitoring the variation of the system energy with respect to the coordinates of the atomic nuclei allows the investigation of hypothetical reaction mechanisms, because one can construct potential energy profiles as a function of proposed deformation pathways that connect reactant states to reaction intermediates (including transition states) to product states. The computation of the geometric, thermodynamic, and spectroscopic properties of these states can also be used to interpret the results of in situ measurements and provide cross-validation between the virtual and experimental approaches. A list of possible outputs that can be gained by performing a DFT simulation is provided in Table 1.^10-21  An example of electron distribution calculated from DFT can be found in Figure 1,²²  which shows the frontier molecule orbital density distribution (highest energy occupied molecular orbital [HOMO] and lowest energy unoccupied molecular orbital [LUMO]) and electrostatic potential (ESP) for three different inhibitors on Cu surface. In such studies, inhibitor efficiencies can be related to properties of the electronic structure. More details regarding such studies will be illustrated later in the section corrosion inhibitors.

Electronic structure methods, by virtue of their rigorous treatment of the fundamental electronic principles underlying chemical bonding, can be universally applied to all states and classes of matter. At the same time, this universality comes at the cost of being computationally expensive and restricted to relatively small clusters of atoms either in isolation or as a periodic ensemble (thus allowing the treatment of crystalline solids and surfaces, in addition to molecules). For the purpose of rapidly evaluating reaction profiles or screening large atomic ensembles, much of the detailed electronic structure information is unnecessary and only the relation between the configuration of atom centers and the total potential energy or its derivative with respect to the atomic positions (the forces on the atomic centers that are utilized to perform molecular dynamics) is required. In these cases, an interatomic potential often provides a better solution if a suitable parameterization can be found. Interatomic potentials use a set of parameterized equations to quickly evaluate the energy of an atom, or atom pairs and triplets, within a given locality based on nearest-neighbor distances, chemical notions of bond order, bond angles, torsions, etc. without resorting to determining all of the details of electronic structure. For systems of interest to corrosion scientists and engineers, however, reliable methods for consistently generating accurate interatomic potentials have not yet been developed, leaving one with the task of devising idealized yet still representative model systems amenable to a full electronic structure calculation. In any case, electronic structure methods are frequently required to construct the reference databases that are then used to parameterize and validate forms and parameters of interatomic potentials.²³ 

DFT is presently the most widely applied first principles method for evaluating reaction mechanisms on materials surfaces, as well as evaluating details of the electronic structure of bulk materials, materials surfaces, and point defects.^24-25  It is also in widespread use for studying the properties of molecular systems. The method was derived from earlier band structure methods applied in solid state physics, and it became adopted in chemical systems when it was noted that it provided a way to go beyond the accuracy of the standard Hartree-Fock molecular orbital approach in quantum chemistry.²⁶  The Hartree-Fock method expresses the electronic structure of a molecule in the mathematical form of a determinant consisting of the individual molecular orbitals, as a means to fulfill the Pauli exclusion principle (electron exchange). It turns out, however, that a single determinant is not sufficient, in the general case, but multiple kinds of determinants can be generated, and the mixing between them lowers the total quantum mechanical energy of the system, a mixing effect known as electron correlation.²⁶  Solving for multiple determinants is too computationally arduous for systems of significant atomic complexity, and it was discovered that DFT provides a means to obtain an equivalent level of accuracy without requiring the same magnitude of computational effort. This pairing of the physical and chemical theoretical communities generated new innovations at the heart of DFT. DFT solves the one-electron Schrodinger equation (the Kohn-Sham equations), which is written as: where the first term is kinetic energy, is potential due to nuclei, is Hartree potential due to the effective potential of all the electrons, and  is exchange-correlation potential.

Innovation in density functional theory is ongoing because, despite its widespread adoption, the electronic structures and potential energies calculated via DFT come along with a measure of uncertainty that is hard to systematically control.^25,27  This shortcoming is a symptom of the fact that the exact solution for the electronic structure is dependent upon knowing the exact form for the relation between the electron density distribution and the energy contribution that is associated with electron exchange and correlation.²⁸  Various approximations to this exchange and correlation functional are available, having been parameterized to reproduce the properties of certain atomic, molecular or materials systems. Moreover, there are also other functionals that aim to satisfy ab initio conditions, like strongly constrained and appropriately normed (SCAN) semilocal density functional,²⁹  as well as higher order hybrid functionals.³⁰  However, while these functionals do make differences, even when benchmarking data are available, the uncertainty when the method is applied to any new molecular or materials system cannot be systematically deduced.³¹  For this reason DFT studies should be applied in combination with experimental characterization techniques (i.e., ICME) so that extensive validation between the methods can be helpful in clarifying the limits of applicability, and the results of DFT serve as interpretive tools or for rapid screening.^32-33  As reviewed in the next section, the exchange-correlation functional is by no means likely to be the most significant cause of error in the DFT simulation of complex systems where corrosion happens. At the same time, it is essential that authors who publish studies that use DFT modeling include in their papers all of the details regarding how the calculations are performed, including choice of functional, the k-point sampling, energy cut-offs and/or basis sets, pseudopotentials, etc. Likewise, the final optimized geometries of the atomistic representations of the system of interest should be made available through supporting information or some other online, readily-accessible portal.

As with any scientific investigation—experimental or theoretical—the decision to use DFT should be made with an awareness of the risks and benefits associated with the technique. One of the key benefits of DFT is its ability to characterize, with atomic resolution, chemical reaction mechanisms, using a bottom-up approach in which complexity is incrementally introduced. This approach contrasts with most experimental methods that use a top-down method by which attempts must be made to incrementally strip complexity away (for example, consider the effort that has to be made to study the behavior of very pure metal systems and idealized surfaces in pure water environments). This complementarity is a strong argument for incorporating DFT into an ICME approach to corrosion science.^34-36  At the same time, risks can emerge due to the need for abstraction of the physical system into a virtual model that can easily be evaluated via the DFT algorithm, and the hard to quantify uncertainties inherent from the combination of abstraction and the underlying physics associated with exchange-correlation.³⁷  A list of risks and benefits of incorporating DFT into ICME is provided in Table 2.

As the primary outputs of DFT calculations are basic physical quantities, such as a list of equilibrium atomic positions, forces per atom, stress tensors on the unit cell, electron densities, electronic band structure or the electronic potential energy, it is usually necessary to apply some “postprocessing” that connects these basic quantities to their observable and experimentally measurable counterparts. For example, potential energies for different phases of a substance can be converted into relative Gibbs free energies through the inclusion of entropic and enthalpic terms derived from the statistical mechanics treatment of the partition factor, that includes, for example, energies associated with the normal vibrational modes of the material. Similarly, terminal solubilities can be extracted by calculating the energy associated with placing an impurity in a material and, likewise, performing a statistical mechanics analysis.¹¹  As a more complicated example, consider the estimation of overpotentials associated with electrochemical reactions on electrode surfaces by comparing the relative adsorption energies of intermediates associated with reactions like oxygen reduction,³⁸  or the prediction of diffusion constants and chemical rate constants from the determination of transition states and their energies relative to isolated reactant states.^39-40  Figure 2(c) shows the calculated overall activation energy of initial oxygen reduction as a function of electrode potential by combining two elementary steps shown in Figure 2(a).³⁸  The transformation between chemical and thermal modes of energy is intrinsic to almost all processes in materials science, including corrosion. The ability of DFT to provide a first principles calculation of energies associated with different chemical or materials states, therefore, creates numerous possibilities for predicting a diversity of phenomena related to these topics.

DFT simulations require the specification of the position of each atom within the abstracted version of the system of interest, followed by the determination of an electron distribution for each electron associated with the atom in the abstraction (often only the valence electrons). Given that a droplet of water 10 nm in diameter would roughly consist of 17,500 water molecules and 140,000 valence electrons, or that a 5-nm-wide cubical grain of nickel might contain 11,500 atoms and 115,000 valence electrons, there are a large number of degrees of freedom one has to specify for even a very spatially modest attempt to simulate a problem in corrosion science. Furthermore, one must consider that the number of computational operations scales with the number of valence electrons cubed, i.e., order(N³). As a result, DFT calculations are rarely performed for more than a few hundred atoms, and usually far less. For this reason, the scientist examining physical problems using DFT adopts one of two primary approximations: the cluster approximation, or the assumption of periodic boundary conditions.⁴¹  Figure 3 shows the structure configuration when calculating O adsorption on Cu with cluster approximation (Figure 3[a]) and periodic boundary condition (Figure 3[b]), respectively. Either approximation of the practical system of interest will have its own advantages and disadvantages. Close collaboration with experimentalists and characterization experts is an essential component to developing models that maintain a reasonable physical fidelity.

The cluster approximation involves the simulation of a modest, finite number of atoms embedded either within the vacuum or some medium that has certain dielectric properties designed to emulate a solvent or continuum environment. Solids can also be studied in this way by creating a virtual “chunk” of some crystal that includes the defect or impurity of interest, and then applying boundary conditions to the atoms on the edge to constrain the behavior so that those atoms mimic the bulk,^42-44  which is also called embedded cluster model.^45-47 

In contrast, the periodic boundary condition model effectively makes every simulation scenario into its own crystal structure. This method is advantageous for including long-range effects for periodic systems like metals, in which truncation to a cluster would lead to considerable “edge effects” in the atomic structure. This method is widely used for surfaces of solids (in which the models are called “slab models”⁴⁸ ) and also for performing molecular dynamics simulations in liquids and solids.⁴⁹  It should be noted for both of these approximations, if the system is too small, and the energetics and even band structures could be wrong.

Both cluster and slab models can be used to simulate ensembles of atoms that can mimic the states relevant to corrosion reactions. For instance, a series of atoms arranged to represent a metal surface using the slab or cluster models can then have various water or other molecules distributed over it to simulate the interaction between solvent and solutes over surfaces.^50-51  Comparing the energy of different arrangements could be a method of determining the likely surface coverage at different values for the chemical potential.⁵²  Using statistical mechanics, one can also include the effects of varying the temperature, pH, and other environmental variables, thus overcoming the once common misconception that DFT calculations are only valid at absolute zero temperature and ultra-high vacuum conditions. Pourbaix diagrams can be theoretically generated from entirely ab initio principles in this way.⁵³  Arrangements of atoms postulated to be metastable states can also be compared to the lower energy states to get a measure of how likely these states will be accessed during a proposed corrosion mechanism, for example, using Boltzmann statistics.¹¹ 

As another example, defect states in materials, such as grain-boundaries or dissolved species like hydrogen or carbon, can be simulated starting with an ideal crystal and then generating the defects by geometric transformations. Using these methods, the energies of different grain boundaries and surfaces can be compared, as well as the energy required to inject an atom, like hydrogen, into the bulk crystal interstices or into the grain boundary. This method allows the quantitative evaluation of mechanisms that require diffusion into the solid state and the interaction of impurities with geometric defects.^54-55  Using various comparative methods and constrained calculations, one can also separate the role of chemical versus mechanical effects in controlling the overall interaction energies.¹⁹  An example is the use of DFT to evaluate the change in zirconium grain boundary cohesion energies when various chemical species are present.⁵⁶  Figure 4 shows that the preferred locations of iodine are Zr-free surfaces and substitutional at a Zr grain boundary rather than bulk Zr. And Figure 5 shows the effect of different impurity elements on grain boundary strength, which shows that other than Nb, Cr, and Fe, all other elements tend to weaken the strength of grain boundary.

Once an atomistic scenario has been constructed to resemble the physical problem of interest using either the cluster method or the periodic boundary systems, the type of DFT calculation must then be selected. Most commonly, a geometry optimization will be performed: the coordinates of the atoms specified in the initial guess—the geometry constructed by the scientist as a “first pass” at what they think the structure of the solid, interface, or molecule should generally look like—are adjusted so as to find a local minimum in the potential energy. By this process, bond-lengths will expand or contract, bond angles and torsions will flex, and/or crystal symmetries and lattice parameters adjust, in such a way as to find the lowest energy structure in the general neighborhood thus indicating the nearest “stable state” for the atomic configuration relative to that initial guess. With information about the various locally stable states, integrating those states and knowledge about the relevant environment into a statistical mechanical model allows understanding of those structures relevant to experiment. The process of geometry optimization can also be used to search for minimum energy pathways, transition state structures, and metastable states with constrained optimization; therefore, this method is often used to elucidate the potential energy associated with given reaction coordinates, and hence to compare and evaluate proposed reaction mechanisms. Whereas the final outputs of the geometry optimization provide physically meaningful states, the states used as initial guesses and the pathway of the optimization method itself are not physically meaningful, in as much as they are by-products of the algorithm used to perform the geometry optimization and the guesswork of the scientist building the model.

As geometry optimization techniques are local in nature, techniques have also been developed to go beyond the limitations of having to only find states close to the initial guess. These methods aim for the goal of searching for the global minimum energy state and involve searching larger swathes of the conformation space of the atoms, molecules, and materials involved in the model system to find the lowest energy states. Metropolis Monte Carlo is an example of this method, as is simulated annealing.⁵⁷  These methods can be combined with machine learning approaches, such as genetic algorithms, to speed up the search.⁵⁸ 

In cases where a reaction pathway is unknown, or one wishes to sample many different conformations for a fluctuating system (such as a solvent, or a material subject to radiation damage) then molecular dynamics can be used to propagate some initial state in time. The molecular dynamics method uses DFT or interatomic potentials to compute the force on the atoms in a given state, and then uses classical molecular dynamics and time-integration to compute the positions of the atoms at a future time step, and then iterates this process for a fixed number of steps specified by the user. It is therefore useful for simulating the time-evolution of a system, and to explore its equilibrium behavior over a statistically well-sampled time period. Another feature that can be coupled with molecular dynamics is the inclusion of constraints that force a certain reaction step to occur (such as a breaking of a bond or the application of a strain rate).⁵⁹  In this case, the variation of the system energy, and the modes for relaxation and the evolution of a response to the constraint event, can be monitored over time. Obtaining adequate samplings for molecular dynamics simulations is a challenging endeavor, however, especially when one is constrained to using DFT only and reliable interatomic potentials are unavailable.⁶⁰  Figure 6 shows the initial and final structure of dissolution of Cu adatom to solvent.⁴³  Here, Taylor tried to explore the dissolution path of metal atom to solvent with DFT. A Cu adatom was moved along z axis (perpendicular to surface) with a step size of 0.02 nm, and at each step, an ab initio MD calculation was first implemented to relax water molecules then a DFT calculation was done to calculate the energy of the system. Figure 7 shows the final results of system energy as a function of bond-breaking distance under different conditions.⁴³ 

In this section, some specific examples will be provided of where DFT has been applied to investigate atomistic and molecular phenomena that have a bearing upon outstanding issues in corrosion science and engineering. This will begin by reviewing the topic of corrosion inhibitor modeling, then move on to the topics of metal oxidation and passivation, localized corrosion, and the modeling of dissolution of metallic materials. In the section following, a forward-looking perspective will be provided on what exciting opportunities exist in the near future for the modeling of corrosion processes with DFT.

Due to the relative simplicity of modeling molecules in isolation using DFT and related electronic structure methods, the quantitative structure activity relationship method has been applied for some time to investigate the properties of molecules found to inhibit corrosion, in order to look for underlying properties that are intrinsic to good inhibitors.^61-66  Obot, et al.,⁶⁷  has given a comprehensive review of the application of DFT on corrosion inhibitors recently. Here there will just be a brief overview. Some of the properties of inhibitor molecules that can be extracted from electronic structure calculations include:

• •

  Molecular orbital structure and energy levels. The frontier molecular orbitals (called HOMO and LUMO) are the molecular orbitals that are most reactive and can donate or accept electrons, respectively, through overlap interaction with, for example, the dangling bonds that appear on metal or oxide surfaces. On one hand, a high value of E_HOMO is likely a tendency of the molecule to donate electrons to appropriate acceptor molecule of low empty molecular orbital energy. On the other hand, a low value of the E_LUMO indicates the tendency of the acceptance of electrons from the d-orbital of the metal. Therefore, the energy difference (E_LUMO – E_HOMO) between these orbitals controls how easily the frontier orbitals can “hybridize” to adapt new bonding relationships with metal surfaces, and hence can also direct how strongly the molecule will bind to the metal surface.^68-69  Figure 1 shows the structure with LUMO and HOMO of different corrosion inhibitors calculated from DFT.²²  In this study, three different inhibitors are calculated with DFT to study their inhibiting efficiency. The HOMO-LUMO gaps for amino pyrazolo-pyrimidine (APPH), 4-hydroxy pyrazolo-pyrimidine (HPP), and 4-mercapto pyrazolo-pyrimidine (MPP) are 5.358 eV, 5.404 eV, and 3.489 eV, respectively. While the order of these value is not exactly the same as the inhibitor efficiency, it does align with experimental observation that MPP, which has much lower energy difference compared to the other two, has the highest efficiency.

• •

  ESP map indicating regions of positive, negative, and neutral charge. These maps, often represented as a color-mapped isodensity surface around the space of the molecule, provide a visual way to determine which parts of the molecule will be most likely, from an electrostatic point of view, to engage with the surface and align in the electrochemical double layer. ESP maps of different corrosion inhibitors are also shown in Figure 1 as the last column.²²  Here, the red represents electron-rich region while blue represents low electron density region. As can be seen from this figure, MPP contains more electron-rich area than HPP, indicating MPP molecules can more easily share or transfer free electron pairs to copper to form firm covalent bonds than HPP. Also, APPH has more blue area that the other two. Therefore, the order of sharing free electron follows MPP > HPP > APPH, which is consistent with the order of the efficiency of these inhibitors.

• •

  Effective molecular volume. This parameter provides an estimate of the site-blocking capabilities of the molecule at the metal/electrolyte interface. The effective shape of the molecule will also determine how effectively the molecules can align in a self-assembled monolayer to form a well-packed barrier to corrosion, although shape is harder to quantify than volume alone.

• •

  Atomic charges and/or partial charges. Various methods are available to estimate the charge on a given atomic center by partitioning the electron density into atom-by-atom quantities. Despite decades of effort and frequent requests by experimental collaborators to provide a value for the atomic charge at a given site, no method has found universal acceptance due to the fact that there is no fundamental definition for such quantity. Partial charges can provide similar information to the ESP map.

• •

  Partition coefficient (log P). The partition coefficient indicates the extent to which a molecular species will pass from one solvent phase into another,¹⁶  and can be extracted from DFT calculation using free energy calculations and implicit solvation models.⁷⁰ 

• •

  Acid dissociation constant (pK_a). The pH sensitivity of corrosion inhibitors is related to both the mode of corrosion at the surface, but also the speciation of the inhibitor molecule as a function of the hydrogen ion activity. When the hydrogen chemical potential is high (low pH), hydrogen ions can bind to electronegative centers, such as secondary or tertiary amines, or, vice versa, hydrogen ions can be released at low chemical potentials to form basic analogs of the inhibitor. The free energy for this reaction can be calculated using various techniques and the pK_a estimated. An example is the isodesmic method which computes the free energy of proton exchange between the molecule of interest and a reference acid, for which the pK_a is already known.¹⁸ 

• •

  Fukui indices. A relatively new feature that has been used to analyze inhibitor performance is the Fukui index. Fukui indices have been described as follows: “Fukui indices are, in short, reactivity indices; they give us information about which atoms in a molecule have a larger tendency to either lose or accept an electron.”⁷¹  The concept is therefore analogous to the work function of a metal, but applies on an “atom by atom” basis for the entire molecule. Hence, the Fukui indices can be used to determine the relative strength of different potential attachment sites in a molecule and to compare different molecular derivatives to one another.

Although a great deal of molecular information can be obtained by performing an electronic structure calculation on the molecule in isolation, in order to understand its interaction within the context of corrosion, it is essential to study the molecular-level interactions with oxidized, partially oxidized, or bare metal surfaces and their defect states. Kokalj has warned against relying on molecular properties alone to infer corrosion protection ability.⁷²  Figure 8 shows an example of DFT calculated optimized structure of different molecules on Cu(111) surface.⁷³  However, calculations including inhibitor molecules and a surface are more complex, as there are several possible configurations for the arrangement of inhibitor molecules on metal surfaces, and varying degrees of oxidation and defect states possible. Progress has been made in this direction,^72,74-81  although it is unlikely that the same degree of throughput will be possible due to the computational effort required to solve the electronic structure problems for molecules interacting with metal surfaces as compared to molecules in isolation. In this regard, the development of reliable interatomic potentials will be critical, and some work of this nature has already been explored by Khaled.^82-84 

As oxides are critical to metal stability under aqueous and atmospheric conditions, the ability to simulate the nature of oxygen adsorption and the formation of multilayers on metal surfaces with various compositions and defect states promises to be exceptionally valuable to providing a first principles rationale and design strategy for corrosion resistance. At this time, progress has been made in terms of simulating atomic oxygen adsorption on low-index metal surfaces, metal/oxide interfaces of particular crystallographic relationships between the oxide symmetry and that of the metal, and performing larger-scale molecular dynamics studies of metal-alloy/oxide interfaces.^85-90 

If one begins by considering atomic oxygen adsorption on a bare metal surface, then the complexity at the atomic level has to be appreciated.⁹¹  Each configuration has its own energy, which will determine how stable the configuration is compared to another. As atoms diffuse around, according to the thermally-driven random walk, the configurations with the lowest energy will be the most stable. DFT can be used to evaluate these potential energies by (a) generating the structure according to some initial guess, (b) performing a geometry optimization to relax the bond-lengths and atom positions according to the quantum mechanical interactions between the atoms and their electrons, and (c) extracting the potential energy from the optimized structure. In the case of partial or complete monolayers, the results can sometimes be compared to ultra-high vacuum surface measurements using x-ray spectroscopy, electron microscopy, low energy electron diffraction or electron energy loss spectroscopy. The interaction of atomic oxygen with metal alloyed surfaces has seen a lot of attention from the DFT community due to its relevance for oxygen reduction reaction electrocatalysis and other surface reactions in heterogeneous catalysis, driven by the intense effort applied over the past decade or so to improve upon precious metal catalysts and electrocatalysts for application in the chemical processing industry as well as fuel cell systems.²⁴ 

Not only are there different configurations for the atomic oxygen on the surface, but there are different coverage levels that should be explored, and configurations within those coverage levels.^92-93  Thus, the atoms may arrange differently at very dilute coverages, to when the coverage is increased to 0.25 ML (a quarter monolayer), 0.5 ML, and 1.0 ML. Beyond 1.0 ML, ultra-thin oxide layers will be formed, and then, beyond that, the bulk oxide will begin to be expressed. Figure 9 shows the difference of O adsorption energy with different oxygen coverage to both Ni(111) and Ni-Cr(111) surfaces.⁹²  Depending on the stress imposed on the oxide by the underlying metal crystallography, the oxide film structure may or may not conform to the most thermodynamically stable bulk oxide phase. The possibilities are numerous and require careful evaluation, an evaluation that can be made within DFT. The effect of geometrical defects on the surface can also influence the nature of the oxygen adsorption and arrangement, and so features such as steps, vacancies, and kinks have also been investigated using DFT. To assist in scaling up atomistic DFT calculations, coarse-graining methods like the cluster expansion can be used to take into account adsorbate-adsorbate interactions.⁹⁴ 

One of the current challenges that faces DFT calculation of the metal/oxide interface is that metal-oxides can have electronic states that the current generation of DFT methods finds hard to accurately describe. This is due to the fact that DFT performs well for metallic-like systems where electrons are largely delocalized, but when the metallic orbitals are localized, as occurs during the metal to oxide transition, then the exchange-correlation functional approximations can break down. There are several proposed methods that have been applied to try to overcome this problem, like local spin density approximation and the Hubbard U term (the LSDA+U method),^95-97  meta-GGAs (generalized gradient approximation),⁹⁸  and screened hybrid functional,^99-100  but none have obtained universal acceptance, and they usually make the calculations more computationally complex, and thus time-consuming to solve. An example where methods for incorporating these localized electronic structure corrections includes the work by Olatunji-Ojo on the growth of oxide films on nickel surfaces,¹⁰  and a similar method was applied to investigate the growth of oxides on plutonium and technetium metals.^101-102 

Beyond the study of these ultra-thin oxide films on metal surfaces, there have also been studies performed on the nature of point-defects in bulk oxides, relevant to the point-defect model for oxide film growth, as reviewed by Uberuaga, et al.¹⁰³  As of yet, however, there has been no comprehensive integration of these works to provide a complete first principles model for passivity but it seems that, given a concerted effort, such a program could be accomplished, at first for a specific single-component metal, and then for binaries and more complex alloy compositions. The first principles outputs, such as activation barriers for point-defect migration, and the stress-strain relations for different oxide structures over the metal substrate will need to be integrated into a continuum model for the oxide growth, possibly in combination with a phase-field model that includes mesoscopic parameters,¹⁰⁴  such as the grain-boundary interfacial energies and the contributions of dislocations as pipelines for point-defect transport.¹⁰⁵ 

Geometry optimization often yields unexpected insights: starting with a naive expectation of an interfacial structure, say, assuming epitaxy or a certain surface termination, the potential energy surface may actually indicate that the expected structure was energetically unfavorable, and relaxation will indicate the closest energetically preferred structure (which may not be the global minimum in energy, as the simulation is biased by the choice for the initial guess). Accordingly, geometry optimizations provide an alternative means for characterization and discovery. An example can be found in the work of Olatunji-Ojo, et al., in the investigation of structural relationships for the Ni/NiO interface,¹⁰  which provides a complement to the scanning tunneling microscopy (STM) investigations of Seyeux, et al.¹⁰⁶ 

Likewise, a DFT study of partially oxidized Mg surfaces revealed that the preferred adsorption site for oxygen was subsurface, not exterior to the surface (as is usually expected), and, not only that, but the potential energy analysis indicated that the oxygen atoms on Mg preferred to cluster rather than repelling one another.⁵²  This second behavior was also unexpected, and was further investigated by Todorova, et al.,¹⁰⁷  and they found that there is a critical height of 1.6 Å (0.1 nm) above the surface at which the interaction becomes repulsive, and that the second row electronegative elements adsorbed on Mg(0001) exhibit an unusual behavior of reduction of the work function with increasing surface coverage. Later, Poberžnik and Kokalj¹⁰⁸  studied the attractive interaction between O and Al and found that adsorption induced decrease of the work function is not necessary for the attractive interaction to occur. They found that the attractive lateral interaction between negatively charged adatoms on Al surfaces are a consequence of an electrostatic stabilization that stems from an interplay between ionic interactions and geometric effects, and the model should be generally applicable provided that (i) the adsorption bonding is sufficiently ionic and (ii) adatoms are sufficiently small to come close enough to the metal surface. These works provide a foundation for further work on the mechanisms of metal oxidation, and, in time, should produce insights as to how stainless lightweight alloys may be developed from first principles.

A first principles molecular dynamics study was performed by Motta, et al., in 2012,¹⁰⁹  and used to demonstrate the interplay between hydrogen bonding, proton transfer (Grotthus mechanism), and defect states on AlOOH (boehmite) surfaces exposed to water. The results demonstrated reactivity of defective surfaces with an explicit solvent, and provided a basis for future work examining interactions between aggressive ions such as H⁺, OH⁻, and Cl⁻ interacting with aluminum metal and its alloys.

Localized corrosion remains one of the most challenging fields of corrosion to model. Various causes have been attributed to localized corrosion, such as the galvanic acceleration of corrosion due to intermetallic particles, the role of line and point-defects, embrittlement and localized thinning of the oxide film, the coalescence of vacancies to form voids at the metal/oxide interface, or influence of aggressive ions.^110-115  Macdonald, et al., proposed the Point Defect Model^116-119  about 40 years ago, which suggests that it is the condensation of cation vacancies (flux of cation to solution larger than flux of cation vacancy to metal matrix) at the metal/oxide film interface that leads to breakdown. Frankel, et al., recently gave a short overview of the current status of localized corrosion¹¹³  mainly from experimental prospective and proposed a new model on pitting stability^120-122  considering both oxide breakdown and pit growth stability. With the help of aberration-corrected transmission electron microscope and a fast and precise super x-ray energy-dispersive spectrometer, Zhang, et al.,¹²³  have recently found that nano crystals (NCs) are embedded in amorphous oxide film, and the NC/amorphous interface provides a path for diffusion of Cl⁻, which allows them to penetrate oxide film and reach metal surface to cause breakdown. While progress has been made in the last few decades, detailed mechanisms of passive film breakdown and pit initiation have not been provided yet. DFT has been a useful tool on studying the atomic scale interactions at surface and interface and potential mechanisms causing localized corrosion, including the interaction among different ions,^124-126  between ions and oxide film,^112,124-131  between ions and bare metal surface,^90,132-134  as well as the effect of defects.^{112,131,134-136}  A few examples will be shown below.

Bouzoubaa, et al., in 2009 used a periodic DFT model to explore the potential for Cl⁻ to displace OH⁻ at a hydroxylated NiO(111) terrace surface.¹²⁴  The researchers used geometry optimization to investigate the thermodynamic tendencies associated with adsorption and subsurface insertion of the Cl⁻ ion into the surface. The study was performed as a function of chloride coverage, where it was found that Cl⁻ is only energetically favorable to adsorb at low coverage, whereas subsurface insertion was energetically favorable at higher coverages. The findings of the study, however, were inconclusive and required further thermodynamic and kinetic analysis. In a follow-up study, the authors reexamined the interactions on a stepped surface, the hydroxylated NiO(533) surface. Fragments of isolated Ni(OH)Cl and Ni(Cl)₂ clusters were formed upon geometry optimization,¹¹²  as shown in Figure 10, indicating that the Cl⁻ may play a role in the dissolution of the oxide films. However, again a more detailed thermodynamic or surface science model (such as deriving a kinetic equation or adsorption isotherm) would be required to embed the DFT values into a localized corrosion context.^112,124-125 

In an interesting paper from 2014, Schmidt, et al., presented a series of DFT calculations of metal/interfacial systems comprised of metal slabs with explicit water molecules, ions, and alloying elements.^134,137  A “progression of geometry optimizations” was shown illustrating differential reactivities of cupric and ferric chloride on Ni(111) with defects (steps)—although, as stated in the above section on DFT methods, the progression of steps in a geometry optimization should not be interpreted as a step-wise reaction sequence, but rather as an algorithmic refinement to find an optimum energy state using steepest descents or conjugate gradients. In this way, the DFT has been used in more of a notional sense to convey mechanistic and energetic preferences, rather than absolute fundamental information regarding hypothetical reaction mechanisms.

In 2015, Liu, et al., explored the influence of chloride on potential pitting surface chemistry of aluminum metal.¹³³  The adsorption of Cl⁻, O atom, and O₂ were computed, and then the potential for Cl⁻ to disrupt oxygen monolayers was calculated. Figure 11 shows the optimization structure of surface with different Cl^– coverage. It can be seen from Figures 11(c) through (e) that Cl⁻ pulled up Al atoms below them, causing considerable disruption to the first Al layer. The average height difference within the first Al layer increases from 1.87 Å to 2.28 Å (0.187 nm to 0.228 nm) with the coverage of Cl⁻ increasing from ¼ ML to ¾ ML. However, when Cl^– coverage is 1 ML (not shown in figure), ions move upward about 1.9 Å (0.19 nm) without affecting the oxygen layer. It suggests that the potential for Cl⁻ to destroy the oxygen monolayer was weak. The formation of AlCl₃ and Al₂Cl₅ species and substructures was shown to be possible, however, as shown in Figure 12. Extension of this work in future efforts, such as to compute the surface coverage as a function of Cl⁻, pH, and temperature using first principles thermodynamics could be very illuminating.

Islam, et al., in 2009 published an examination of the condensation of atomic vacancies at the oxide/alloy interface between alumina and TiAl(111).¹³⁶  The formation energies and relative formation energies for vacancies at the interface were calculated and compared to the bulk vacancy formation energies. The mechanism that was revealed in this paper involved the formation of 2D clusters of vacancies.

Besides what has been described above, the use of reactive force field potentials has also been applied to simulate localized corrosion phenomena, such as the chemical breakdown of CuO films on Cu by chloride,¹³⁸  and the influences of interfacial electric fields upon the initial oxidation of Ni(111) in water. Although these are not within the scope of DFT, it should be noted that DFT is frequently used to construct the potentials used in such simulations. As discussed in the introduction to theoretical methods, force-field based molecular dynamics simulations can go beyond DFT in exploring numerous types of stable sites that can form at ambient or elevated temperatures, thus going beyond the limits of the time-consuming and guess-dependent geometry optimizations performed in a DFT calculation.

To design corrosion-resistant alloys (CRAs), it is essential to understand the detailed mechanism of dissolution, which is part of the process of localized corrosion. Dissolution of metallic materials is usually controlled by multiple factors including the surface configuration, crystallographic orientation, alloying elements, solvent environment, and applied potential. With the development of computer technology, DFT has emerged as a promising method to study the atomic scale mechanisms of dissolution of metallic materials in the last two decades. These models not only allow the study of fundamental mechanisms of dissolution, but also can allow one to study a wide range of materials to design CRAs that otherwise would not be possible for experiments. In this section, the application of DFT on dissolution of metallic materials will be reviewed.

Qualitatively, DFT can be applied in two ways. (1) Study the corrosion potential of alloys and find the most promising elements to increase the corrosion resistance of alloys. According to Fermi-Dirac function, a lower work function implies a lower corrosion potential. Han, et al.,^139-140  calculated the work function associated with water adsorption on γ-Fe(111) surface with 10 different alloy elements, and found that most of the alloying elements—especially Cr, Mo, and Si, but not Ti and Nb—can mitigate dissolution corrosion of the γ-Fe substrate in an aqueous environment. They also found that when Cr combines with interstitial atoms, especially in 2Cr-N and 4Cr-C co-doped surfaces, it leads to a weakening of the surface corrosion resistance. Greeley and Nørskov^141-142  proposed a simple thermodynamic formalism for estimating the dissolution potential of solute metal atoms in various host metals referenced to the dissolution potential of the solute in its pure metallic form using DFT. The model assumes a grand canonical ensemble with saturated host ion in the liquid, thus the energy of the system can be expressed as a function of the electrode potential difference (U-U₀) between the alloy and pure bulk solute element. Therefore, the energy change of the system before and after dissolving a solute atom can be calculated, which is a function of U-U₀ as well. With this model, they studied the dissolution of various solute/host combination binary alloys, with solute embedded in the surface,¹⁴¹  as isolated adatoms, dimers, and kink structures.¹⁴²  Figure 13(a) shows the structures used in these calculations, including isolated adatoms (a1), dimers (a2), and kink structures (a3), and Figure 13(b) shows the calculation results for each kind of structure. It can be seen from Figure 13(b) that the deposition/dissolution potential shifts are, on average, very negative over the entire set of admetal/substrate pairs for isolated atoms and dimers, although there is an average shift to higher deposition/dissolution potentials for dimers compared to isolated atoms. However, for kinked structures, there is a higher shift of deposition/dissolution potentials when adatom is noble metal. These results can help guide the experimental searching of alloying element. This model has later been successfully applied to other systems,^142-147  including Al,¹⁴³  Co,¹⁴⁴  Ni,¹⁴⁵  Pt,¹⁴⁶  and Mg.¹⁴⁷  (2) Study the mechanism of dissolution, e.g., figure out each elementary step during dissolution process. There are only a limited number of research studies^148-150  that have been done in the area of metallic material corrosion. But it has been successfully applied in other areas, like dissolution of minerals^151-154  and salts.^21,155-156  Liu, et al.,¹³³  studied the effect of Cl⁻ coverage on dissolution of Al, and found a coverage between 7/9 ML and 8/9 ML is the critical coverage of Cl⁻ to the active dissolution of Al in pitting. They found that the two AlCl₃ substructures formed in the 7/9 ML case are pyramids with similar bonding lengths and angles, whereas the AlCl₄ substructure formed in the 8/9 ML case is a tetrahedron with its four sides as well as its six angles nearly equal to each other. Connecting the details of these geometries and configurations to mechanisms of corrosion and predictive performance has yet to be achieved.

Quantitively, DFT can also be used to provide fundamental parameters to higher scale modeling, like kinetic Monte Carlo and continuum scale modeling. Ma, et al.,¹⁵⁷  recently published a paper on obtaining anodic dissolution curves based on DFT calculations of work function and surface energy density, which builds a bridge between atomistic modeling with experiments. Figure 14 shows the comparison of their calculation results and experimental data.^158-159 

The dissolution rate of a metal can be obtained via Tafel-type equation:¹⁶⁰ where n is the number of electrons involved in the electrochemical reaction; k is a reaction constant; and F, R, and T are Faraday constant, gas constant, and absolute temperature, respectively. ΔG_ch is the chemical activation energy for a metallic ion to escape from the metal lattice and dissolve into the solution, and ΔG_ele is the modification of dissolution activation due to applied potential and double-layer potential.

An accurate calculation of the chemical dissolution activation energy ΔG_ch is the key to get quantitative results of dissolution rate of metallic materials. However, because it is difficult to determine the transition state of dissolution from metal atom on surface to ion in solution, the activation energy has been approximated or determined from the following ways with the help of DFT:

• •

  Bond energy. Dissolution activation energy can be estimated from bond energies as a dissolved atom needs to break the bonds with its neighbor atoms before dissolving into solution. Dissolution activation energy can be estimated as ΔG_ch=∑_i=1ⁿ E_i, in which n is the number of neighbor atoms of the dissolved atom and E_i is the bonding energy with each neighbor atom. Researchers have either set the bond energy to a constant number^161-162  or calculated the bond energy between different elements^163-166  with DFT or interatomic potentials.

• •

  Surface energy. Dissolution activation energy can also be estimated from surface energy as ΔG_ch = βγ,^147,167  where γ is the surface energy and β is a coefficient. According to this approach, Song, et al.,¹⁶⁷  showed that the dissolution rate of the (⁠⁠) plane is 20 times faster than that of (0001) plane for pure Mg. DFT can be very useful to calculate surface energy of alloys. A lot of papers^{88,143,147,157,168-170}  have been published recently regarding the calculation of surface energy of alloys using DFT with the effect of alloy composition and crystallographic orientation. For example, Deshmukh and Sankaranarayanan⁸⁵  studied the effect of impurity and crystallographic orientation on surface energy of Mg alloy with DFT, and found that surface (0001) has a lower energy compared to surface (⁠⁠) and (⁠⁠), and that the existence of Al, Zn, and Y impurity elements on the Mg surface reduces the surface energy, while adding Ca increases the surface energy.

• •

  Direct calculation with constrained minimization approach. Taylor¹⁷¹  used DFT to investigate the dissolution of Cu by lifting a Cu atom from a nano particle and fixing the position of the lifted atom along the z axis. He also studied the effect of applied potential and the impact of binding water molecules on the dissolution activation energy, as shown in Figures 6 and 7. Pham, et al.,¹⁴⁸  studied the dissolution of Fe adatom from Fe (110) surface slab and the effect of Mo and Cr by elevating the adatom along the z axis, in vacuum and also Fe(OH)₂, which mimics the environment of water. Similarly, Wang, et al.,¹⁵⁰  studied the Mn ion dissolution from MnS with DFT. The transition states of each elementary reactions of dissolution of Mn ion were found by doing constrained geometry optimization first along the y axis and then along the z axis.

While the methods above can provide some insights on the dissolution of metallic materials with the help of DFT that otherwise would not be possible with experiments or other modeling methods, there are still some limitations. For example, although it is reasonable to estimate the dissolution activation energy from either bond energy or surface energy, these methods are more suitable to calculate the relative dissolution rate either for different impurity elements or crystallographic orientation rather than the absolute dissolution rate. Taylor, et al.,^148,150,171  made one step forward to calculate the activation energy with constrained minimization approach; however, the detailed mechanism during a dissolution process is still not fully understood and it is still a hard problem to do an exhaustive search and determine the identity of the ion in the transition state. Recently, dissolution has been studied by methods like thermodynamic integration (blue moon–ensemble) and metadynamics¹⁷²  in other fields, like atmospheric chemistry (salts),^21,155  and batteries.¹⁵⁶  With the development of computation models as well as computing speed, more insights can be obtained on the dissolution mechanism of metallic materials with DFT, which can lay the foundation to design CRAs.

Besides what has been described above, DFT can also be useful in other applications of relevance to corrosion science and engineering:

• I.

  Thermodynamic assessments. Thermodynamics of a system provides a foundation for all of the other properties of the system. CALPHAD is a widely used method to evaluate the thermodynamic properties of systems. It allows the extrapolation to higher order multicomponent systems by combination of lower order systems and has been successfully applied to materials in different areas. It even allows extrapolating to systems far away from the conditions (composition and temperature) that experimental data were obtained. However, there are certain conditions it is difficult, if not impossible, to obtain experimental data for a system. For example, there is no Cu with body-centered cubic (bcc) structure that exists in nature, but it is needed for the study of early precipitation of Cu in bcc Fe.¹⁷³  In such circumstances, DFT has emerged as a useful tool to complement experimental data and therefore provides thermodynamic parameters needed to build thermodynamic models otherwise not possible. Kormann, et al.,¹⁷⁴  detailed the limitations of CALPHAD models based purely on experimental data and how DFT calculations can supplement. It has been successfully applied in a couple of areas, like nuclear materials^175-178  and thermoelectrics.^179-180  Figure 15¹⁸⁰  shows an example of calculated mixing of enthalpy of (Cr, Mn)₃Si structure by DFT. While not related to corrosion-resistant materials, this serves as an example showing how these calculations can help build CALPHAD models. For corrosion-resistant materials, Liu, et al.,^181-182  built the CALPHAD models for Al-Co-Cr-Ni quaternary system and all of the subsystems by combining DFT results and experimental data. This alloy is used in the protection of high-temperature Ni- and Co-based superalloys operating in corrosive environments such as gas turbine hot sections. More applications of DFT-based CALPHAD are foreseen in designing corrosion-resistant materials.

• II.

  DFT also plays an important role in multiscale modeling (MSM). Currently, there is still a gap between DFT calculations and experimentally measured properties. This gap can be filled by machine learning techniques. However, while powerful, machine learning requires a vast amount of experimental data for a complex system, and usually only works well within the training experiment conditions but may not be reliable for predicting properties under new conditions. A multiscale model that is built bottom up from first principles to continuum scale would be ideal, but no such model currently exists. Multiscale modeling of localized corrosion has been recently reviewed,¹¹⁵  in which the authors detailed the current state of multiscale modeling for the various processes associated with a localized corrosion event, such as formation and breakdown of oxide films, the role of aggressive ions, metastable pitting, and pit or crevice propagation, etc., as well as the proposed improvements that can be made toward multiscale modeling. Taylor, et al.,³²  reviewed the development of ICME in corrosion-resistant materials. Cole and Hughes¹⁸³  reviewed the multiscale modeling of coating and protection systems of CRAs. MSM will benefit from the development of more high-resolution characterization experimental techniques, which will show more lower length scale experimental observation, and allow the validation of atomistic scale models.

• III.

  Another application DFT can be useful in corrosion-resistant materials is the prediction of structures, especially for the surface or interface between material and solvent environment. Finding the most stable structure of a large assembly of atoms is a very challenging problem. Only fully quantum-mechanical calculations are sufficiently reliable to deliver the level of accuracy required. Many global optimization strategies for performing first principles crystal structure prediction have been developed, including ab initio random structure searching,¹⁸⁴  minima hopping,¹⁸⁵  and generic/evolutionary algorithms.^186-187  A few review papers on structure optimization by DFT have been published.^188-190 

• IV.

  DFT can also assist the interpretation of experimental data. This includes both atomic scale experimental data¹⁹¹  by themselves and the underlying mechanism that result in the measured macroscale experimental data. In STM experiments, both oxygen vacancies and hydroxyls on TiO₂ surfaces show up as protrusions on dark rows, which made it difficult to distinguish them. With DFT calculations, Wendt, et al.,¹⁹¹  confirmed that hydroxyls were the brighter ones, which allowed accurate determination of oxygen vacancy concentrations. With DFT, Rashkeev, et al.,¹⁹²  found the important role of hydrogen in pit initiation in aluminum alloys, and Carling and Wahnstrom¹⁹³  found that it is anharmonicity for the lattice vibration rather than the presence of divacancies that explains the non-Arrhenius temperature dependence of the vacancy concentration in Al and a re-interpretation of experimental vacancy data is necessary.

To conclude, DFT has already begun to be applied to illuminate core processes that occur during the corrosion of metals and their alloys. The computational restrictions that enforce certain idealizations and constrain the number of atoms that can be used in a simulation require modelers using DFT to be creative in simulation of corroding systems without introduction of excessive artifacts and loss of practicality. Already the analysis of surfaces, interfaces, and the study of various configurations associated with various elements of corrosion and corrosion control practice have begun to impact sub-disciplines of corrosion science and engineering, such as corrosion inhibition, localized corrosion, and metal dissolution. At this juncture, it is believed that the extent of such studies will expand greatly over the coming few years, providing unprecedented insights into the fundamental mechanisms of corrosion. From this fundamental insight, new and innovative approaches to corrosion prevention and control can be expected to evolve.

The author acknowledges the kind invitation to submit this review article to CORROSION journal from Dr. John Scully; and helpful comments from Drs. Gerald Frankel, Michael Francis, and Narasi Sridhar. This work was supported as part of the Center for Performance and Design of Nuclear Waste Forms and Containers, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0016584.