InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Materials Science » "Adhesives - Applications and Properties", book edited by Anna Rudawska, ISBN 978-953-51-2784-0, Print ISBN 978-953-51-2783-3, Published: November 23, 2016 under CC BY 3.0 license. © The Author(s).

Chapter 14

Ab-Initio Modeling of Adhesive Behaviors at Material Interfaces

By Jun Zhong
DOI: 10.5772/64904

Article top


(a) (Left) An optimized geometry for a fcc-Al/α-Al2O3 interface determined by the UBER calculations; (right) the final relaxed structure: Small spheres are Al atoms, and large spheres are O atoms. The view direction is along [12¯10]. (b) Side-view of electronic structures at an interface: Two slices through the electron localization function (ELF) for an fcc-Al/α-Al2O3 interface taken along (101¯0) and (112¯0) planes, indicating four of hexagonal close-packed oxygen layers in oxide (bottom) and all five atomic layers from an Al slab (top).
Figure 1. (a) (Left) An optimized geometry for a fcc-Al/α-Al2O3 interface determined by the UBER calculations; (right) the final relaxed structure: Small spheres are Al atoms, and large spheres are O atoms. The view direction is along [12¯10]. (b) Side-view of electronic structures at an interface: Two slices through the electron localization function (ELF) for an fcc-Al/α-Al2O3 interface taken along (101¯0) and (112¯0) planes, indicating four of hexagonal close-packed oxygen layers in oxide (bottom) and all five atomic layers from an Al slab (top).
Side-view of geometry at the Al(111)/α-Fe2O3(0001) interface.
Figure 2. Side-view of geometry at the Al(111)/α-Fe2O3(0001) interface.
The final Al(111)/α-Fe2O3(0001) interfacial geometry optimized by the DFT.
Figure 3. The final Al(111)/α-Fe2O3(0001) interfacial geometry optimized by the DFT.
Main sequence of compressing the Al(111)/α-Fe2O3(0001) interface without lubricants.
Figure 4. Main sequence of compressing the Al(111)/α-Fe2O3(0001) interface without lubricants.
The stress perpendicular to the interfacial area as two slabs compress together.
Figure 5. The stress perpendicular to the interfacial area as two slabs compress together.
The number change of Al–O and Fe–O bonds in the R–L region during slab compression.
Figure 6. The number change of Al–O and Fe–O bonds in the R–L region during slab compression.
The main sequence of stretching the Al(111)/α-Fe2O3(0001) slabs.
Figure 7. The main sequence of stretching the Al(111)/α-Fe2O3(0001) slabs.
The number change of Al–O, Fe–O and Al–Al bonds in the new R–L region during the slab tension.
Figure 8. The number change of Al–O, Fe–O and Al–Al bonds in the new R–L region during the slab tension.

Ab-Initio Modeling of Adhesive Behaviors at Material Interfaces

Jun Zhong
Show details


In this chapter, density functional theory (DFT) is employed to identify the essentials of adhesive formation at the Al/Ceramic interfaces. Analyses of electronic structure in the DFT used for atomic bond evolution at interfaces, indicate that, adhesion energy mostly dominated by surface energy, may lead to formation of adhesives, which bind two type bulk-materials strongly together, especially for those polar systems involving cubic oxides, carbides and nitrides. In addition, ab-initio molecular dynamics (AIMD) based upon the DFT is adopted, for example, in the simulation of chemical reaction between two contacting reactive slabs: pure aluminum and iron-oxide. This may provide an insight into the dynamical formation of an adhesive (amorphous Al2O3-texture) occurring between two nascent material surfaces if lubricants are not present or are insufficient. Such a texture may bond onto a hard-roller surface as a protective thin film to resist the elevated temperature.

Keywords: density functional theory, ab-initio molecular dynamics, Al-Ceramic adhesives, amorphous Al2O3-texture, adhesion energy, surface and interface energies, adhesive transfer of materials

1. Introduction

Adhesives have been widely used to improve mechanical properties of materials such as hardening, deformation stability of surfaces or thin films for a variety of technological applications [16]. A significant effort focused on exploring nanomechanisms at interfaces where two solid surfaces contacted. For examples, adhesives were critical to successful manufacture for semiconductor circuits. Each deposited layer must adhere strongly to that beneath it to survive thermal stresses. Similarly, a recent report told that, wear of materials usually resulted in a very high economic loss to every country, for instance, economic losses due to the wear would amount to about 1.0–2.5% of the gross national product (GDP) in the USA [7]. Therefore, adhesives may have a tremendous effect on the US industry because it determines the lifetime of many working parts, from components of automotive engines to airplane brakes. A particular example is the aluminum-alloys processing which is also relevant to many other metallic alloys, that is, the Al-alloy products are manufactured through a series of processing steps that begin with casting an aluminum ingot. In order to reduce ingot thickness, it passes through a series of rolling operations under high loads at speeds up to 4500 ft/min. In this case, there will take place adhesives (the melting aluminum bonds to steel surfaces as chemical additives). On the positive side, many machining operations rely upon adhesive wear processes (from oil drilling to polishing a Si-wafer). The processing optimization is critical to commercial success. While on the negative side, an interest in dry machining is rising because oils used in the machining are likely evaporated due to frictional heating, and such a vapor impact on the long-term health of machine operators is increasingly concerned. Therefore, the better understanding of adhesion mechanisms may result in positive impacts on many industrial manufacturing.

In this chapter, we will introduce our previous studies of adhesives on aluminum and other alloy surfaces: (1) Adhesion is of particular interest. Although a work has reported the adhesion character at an interface between two pure materials (aluminum-diamond) [8], we will discuss an application of electronic structure to this field between aluminum and various ceramic compounds, to determine factors (e.g., an adhesion energy) that may qualitatively control the adhesion strengths at interfaces because this discussion may approach more to the reality even if such computed factors often show too large deviations from some continuum-scale (finite element) calculations. In addition, some particular adhesives, which bind two type bulk-materials strongly at the Al/Ceramic interfaces, are discussed. (2) Another simulation from ab-initio molecular dynamics (AIMD) based upon density functional theory (DFT) is present to bond breaking and formation, and adhesive transfer of materials over an interfacial area where two different material slabs (with pure aluminum and iron oxide) are compressed and then stretched out. In this study, contact of two reactive slabs may lead to a strong exothermic reaction to form a dynamical thin film bonding onto a hard-roller surface to prevent the surface from the melting.

In the following discussions, although we mainly focus on aluminum/aluminum-oxide samples, without loss of generality, we believe what we discuss here is also relevant to many other materials in which adhesion may occur.

2. Methodology

2.1. Density functional theory

Macroscopic materials consist of an enormous huge number (~1023 per cm3) of nuclei and electrons, and these particles interact in a complex fashion. With the development of quantum mechanics in the early 1900s, the law governing those interactions was made clear, that is, the so-called many-particle Schrödinger equation. In order to solve this equation efficiently, the Born-Oppenheimer (BO) approximation [9] was usually invoked to separate the ionic degrees of the freedom from electronic ones. Even so, an exact solution for this many-particle problem was still impossible. In the early 1960s, density functional theory (DFT) was proposed as an alternative to the wavefunction-based description of quantum mechanics. The core concept of the DFT was to describe the state of a system not with the complex all-electron wave function (AE-WF), but with the much simpler charge density. Consequently, the DFT may drastically reduce the degrees of the freedom needed to describe the system, yet in principle maintain an accurate depiction of the physics. In recent years, the DFT has become the method of choice for describing large molecular clusters and solid systems.

DFT was based upon theorems developed by Hohenberg, Kohn and Sham [10, 11], which stated that the ground-state properties of a system were functionals of the charge density n(r) alone. This charge density, which minimizes the energy functional under a given external potential: Vext(r) – the interaction between electrons and nuclei, is the ground-state charge density. So the total energy of a system can be expressed as a functional of the charge density, that is,


According to the Kohn-Sham ansatz [12], the kinetic term T[n(r)]¯ is replaced with a kinetic energy of non-interacting electron gas: EKE[n(r)]; the electron-electron interaction potential Vee[n(r)]¯ is replaced with a classical Coulombic interaction: EH[n(r)] – the Hartree energy, plus an additional non-classical term: EXC[n(r)] – the exchange-correlation (XC) energy. Therefore, Eq. (1) is modified as follows:


The advantage of this modification is that the first two and the last terms in Eq. (2) can be dealt with a traditional manner; while the third term the exchange-correlation energy, whose form is in principle unknown, can be simplified through several approximations. Eqs. (1) and (2) are called the N single-particle Kohn-Sham (KS) equations.

In order to solve the total energy Etot[n(r)] in Eq. (2) more accurately, a set of Lagrange multipliers {εi}, which correspond to the orthonormality of the N single-particle orbitals {φi}, are adopted in an effective Kohn-Sham Hamiltonian operator HeKS, so the ground-state energy Etot[n(r)] can be obtained through minimizing the charge density n(r), see [11, 12]. One advantage of this method is that the number of the N single-particle KS equations increases linearly with the system size and number of electrons.

Even if the Kohn-Sham ansatz is exact in principle for ground-state systems, in practice the approximations are necessary because the exact functional form EXC[n(r)] in Eq. (2) is still unknown. Consequently, an accurate, yet efficient expression for the XC energy is essential. Even though some scientists were devoting themselves to extending the DFT to treat excited states in materials [13], most of current works mainly focus on systems at the ground state.

2.1.1. The exchange-correlation functional

Several approximations for the XC energy functionals have been developed. The most commonly used approximation in the DFT calculations was the so-called local density approximation (LDA) [14], which was based upon the so-called uniform electron gas. The LDA was an exact approximation for a uniform electron gas in a system with a slowly varying charge density [15]. Actually, the charge density n(r) usually behaved like a rapidly varying function of position in some systems. Thus, the corrections to the XC energy based upon the gradient of the charge density were required to reduce the LDA error.

One famous commonly used extension to the LDA was the generalized-gradient approximation (GGA). An example of a GGA functional (PW91) predominantly used in solid-state calculations was that of Perdew and Wang [16]. Another modified expression known as the PBE functional has been suggested by Perdew, Burke and Ernzehof [17]. Comparing the GGA to the LDA, the former has demonstrated to tremendously improve the accuracy of calculations on solids and molecular clusters [15].

2.1.2. Plane waves and supercell

The method of Plane Waves (PW) basis-sets (orbitals), which was more suitable for calculations of periodic solids in the reciprocal space than others on extending its accuracy to any desired convergence, has been widely used for analyzing adhesion properties through calculations of various electronic structures. To simplify, they can be handled by sampling some special k-points in the reduced Brillouin-zone to determine the electronic wave function and the electric potential through the filled electronic bands. Two widely used methods were Chadi-Cohen [18] and Monkhorst-Pack [19] samplings.

In a larger supercell, for the wave vector selection, it is often sufficient to only consider a particular value-point—the Γ point because this point is the original position in k-space, and such a choice would confer significant savings with data storage and computational cost.

2.1.3. Pseudopotential approximation

The Plane Waves (PW) method cannot directly be carried out in the KS theorem since the wave functions of ionic-cores are rapidly varying with many nodes and require enormous PW basis-sets to reach the calculation accuracy. Since chemical reactions are solely related to valence electrons while ionic-core electrons are approximately "frozen" in their original configurations, in practice, to handle the external potential: Vext(r), a scheme of pseudopotential approximation (PP) is provided to just reproduce identical curves of valence wave functions in the place sufficiently far from the ionic-core region (RC). So in the PP framework, the constructed charge densities outside the ionic-core region (RC) must be identical to the true ones, which is called "the norm-conserving condition" [20].

To compromise the computational demanding in practice, a successful concept of "the ultrasoft-pseudopotential (US-PP)" was introduced by Vanderbilt [21]. Blöchl further expanded the US-PP concept by combining ideas from the PP and the LAPW (linearized augmented-plane-wave) methods in a conceptually elegant framework, he called it "the projector augmented-wave method (PAW)" [22]. And then, Kresse established an exact formal relationship between the US-PP and the PAW methods [23].

2.2. Molecular dynamics

Molecular dynamics (MD) is a methodology which depicts motions of a many-particle system using classical Newton’s equations. It has been over 40 years since the first application of MD method created by Alder and Wainwright to a system with hard spheres [24]. The MD method is particularly useful for the study of dynamical properties of materials. More important, this study may help people extract physical insights from the many-particles system.

2.2.1. Equations of motion

In MD simulations, classical trajectories of particles in the real space for a finite size system are generated by solving the Lagrange equations numerically. When choosing the Cartesian coordinates, these equations will become Newton’s equations, that is,


where Vtot(RN) is the total potential energy of a system, mi is the mass of the particle i, Pi is the momentum of the particle i, Fi is the total inter-particle force acting on the particle i, and " " and " ⋅⋅ " are the first and the second orders of the time derivative, respectively. The convenience of solving Eq. (3) should depend on some practical algorithms.

2.2.2. Algorithms

A standard way of solving Eq. (3) is the finite-difference method: Under a given configuration (positions and velocities of all particles) and other dynamical information in a system at the time t, a numerical integration would determine the new configuration of the system at the later time t+Δtt is a time step). Commonly used methods in MD are the Verlet algorithm [25], the Leapfrog algorithm [26] and the Gear predictor-corrector (GPC) algorithm [27]. An ideal algorithm should be simple, run fast and require little memory. It should also permit using a longer time step Δt to hold the whole trajectory of the system as close as possible and preserve conservation laws of momentum and energy. According to these, the Verlet algorithm can be the most widely used one. The leapfrog algorithm is essentially identical to the Verlet algorithm. The GPC method is often more accurate as well as more complicated than those two algorithms.

2.2.3. Statistical ensembles

The widely used statistical ensembles in MD simulations are often the microcanonical ensemble (NVE), the canonical ensemble (NVT), the isothermal-isobaric ensemble (NPT) and the grand canonical ensemble (μVT). Thermodynamic variables in above parenthesis for each ensemble are fixed in relevant ensembles [28]. In MD simulations, the (NVE) ensemble is the best to be realized since all equations of the motion conserve the total energy E and particle’s number N. The constant volume V can be fixed using the free or the periodic boundary conditions; The (NVT) ensemble is very commonly used for practical calculations. Therefore, to control the temperature T in this ensemble, there are a number of ways to simulate it: (i) the simplest one is the velocity scaling method [28], but is crude in MD simulations. (ii) the Nosé-Hoover’s method is a more conceptually fundamental, but has an adjustable parameter of the thermal inertia in practice [29].

For the (NPT) ensemble, Anderson provided an "extended system" method for a constant pressure MD [30]. In this method, a term of "external piston" was added into the Lagrange functional for adjusting a value of external pressure. In another way, Parrinello and Rahman [31] proposed a scheme of "variable cell constant MD" which allowed the simulation box to change the shape as well as the size [31]. This method, so far, is the most general constant pressure method in MD simulations.

The (μVT) is very useful for an open system: In this ensemble, the total particle’s number N will not be conserved, see [32] for details.

2.2.4. Ab-initio molecular dynamics

The MD simulation based upon the force calculations through empirical potential [33] is called the classical molecular dynamics simulation, while that based upon the force calculations using density function theory is called the ab-initio molecular dynamics (AIMD) simulation [34]. Comparing these two, the former can approximately treat a large system (up to tens of millions of atoms) for long time steps (nanoseconds); while the latter can more accurately simulate a small system of hundreds of atoms, and is better at describing reactions involving the change of atomic bonding types.

In principle, the ab-initio molecular dynamics (AIMD) is an extension of the classical molecular dynamics (MD), which includes the first principle-derived potentials for the calculation of the inter-particle forces and the total energy of a system, see [34, 35]. Born-Oppenheimer molecular dynamics

In the AIMD scheme, the total potential energy Vtot(RN) of a system has the same physical meaning as the Kohn-Sham energy (EKS) by means of the Born-Oppenheimer (BO) approximation. Since the Kohn-Sham energy depends only upon nuclear positions and the defined hyper-surface for nuclear motions, the Lagrange functional of the so-called Born-Oppenheimer molecular dynamics (BOMD) has the following formula [35]


where {Λij} are the Lagrange-multipliers. So the minimization of Kohn-Sham energy, EKS, is the constraint to the orthogonal basis-orbitals {φi,j}. And then, Eq. (3) for many-particles can be modified by means of Eq. (4) as follows:




where HeKS denotes an effective Kohn-Sham Hamiltonian operator for a single-electron. Based upon Eqs. (56), a classical molecular dynamics program can be turned into the BOMD program by replacing the energy and the force routines with the corresponding ones from a quantum chemistry methodology [35]. The accuracy of the force calculations used in the BOMD depends linearly on the accuracy of minimizing the Kohn-Sham energy EKS. To realize this, an algorithm of the conjugate gradient (CG) minimization was effectively applied to such calculations, see [34]. Car-Parrinello molecular dynamics

The basic idea of the Car-Parrinello molecular dynamics (CPMD) can be seen as exploring a time-scale separation of "fast-electronic" and "slow-nuclear" motions by transforming it into a dynamical framework of "classical-mechanical and adiabatic energies" [36]. In classical mechanics, forces on nuclei are obtained through the derivative of a Lagrange functional with respect to nuclear positions. The CPMD scheme suggested that, a functional derivative with respect to particle orbitals, which were interpreted as classical fields, yielded the forces on those suitable particle orbitals given by the Lagrange functional. Thus, the CPMD scheme postulates the following Lagrange functional


Like in classical mechanics, nuclear positions, electronic orbitals and the Car-Parrinello’s equations of the motion are modified as the following formulae corresponding to the Newton’s equations [35]




where μ is a "fictitious mass" or "inertia parameter" assigned to the orbital degrees of the freedom. The unit of this "mass parameter" μ is the energy multiplies the time square for reasons of dimensionality. Note that constraints within EKS would lead to "force constraints" in Eq. (8). In practice, these constraints would depend upon both the Kohn-Sham orbitals and the nuclear positions through an overlap matrix of basis-functions—the simulated annealing to search for coefficients of the basis-orbitals that minimize electronic energy. In this terminology, a scheme of "low electronic temperature" or "cold electrons" in Eq. (9) is close to its instantaneous minimum energy: min{φi}EKS[{φi};RN], that is, in this case, it is close to the exact Born-Oppenheimer hyper-surface. Therefore, a ground-state wave function optimized for the initial configuration of nuclei would stay close to its ground state during time evolution if the system is kept at a sufficiently low temperature. A SHAKE algorithm is used to impose constraints on the system to ensure that the KS orbitals remain orthonormal [34]. In practice, it was found that the Car-Parrinello scheme was unnecessary for electronic configurations to be minimized in a coefficient space at each MD time step even though this manner may give errors in force calculations on nuclei. However, it was rather strange to indicate that, the errors of force calculations on nuclei were canceled by the associated errors in depicting electronic motions. According to this, the Car-Parrinello scheme may tolerate force errors on nuclear motions in the AIMD simulation.

In addition, the last term on the right-hand side of Eqs. (5) and (8) would vanish if basis-orbitals {φi} are independent of nuclear positions {RN}. Comparing the BOMD to the CPMD schemes, the former can be viewed as an alternative to the latter, which may give very accurate coefficients of basis-orbitals through minimizing electronic energy at each MD time step.

2.3. An implemented software: VASP

The Vienna Ab-initio Simulation Package (VASP) is a commercial software under the framework of density functional theory (DFT) code that is developed and distributed by a research group at The University of Vienna, Austria [37]. VASP uses a plane wave (PW) basis-sets (orbitals) to expand the single-electron Kohn-Sham wave function and the pseudopotentials for describing electron-ion interactions. Three separate approximations to the exchange-correlation (XC) energy are employed. The local density approximation (LDA) and the generalized-gradient approximation (GGA) are in two flavors: the PW91 functional and the PBE functional [1417]. In addition to these functionals, several additional routines have been implemented to allow for estimation of transitional state energies and geometries (NEB, metadynamics [38]), vibrational properties, analyses of various electronic structure aspects, etc. For more details regarding the VASP code, see [37].

3. Density functional study of adhesion at the metal/ceramic interfaces

A fundamental quantity to influence mechanical properties of a material interface is defined as adhesion energy: Wsep (also commonly referred to as a "work of separation"), breaking interfacial bonds and reversibly separating the interface into two free surfaces, but neglecting diffusion and plastic deformation [3941]. Experimental observations on this dissipative process indicated that, due to plastic deformation, the Wsep was proportional to, but smaller than, the work done to cleave the interface, see [42]. Thus formally, the Wsep is defined in terms of surface and interface energies of respective bulk materials, or to the difference in total energy between an interface system and two isolated surface slabs:


Here σi is the surface energy of slab-i, γ12 is the interface energy, Eitot is the total energy of slab-i, E12tot is the total energy of a system (slab-1 + slab-2) with an interface, and A is the whole interface-area.

3.1. Results and discussions

Available analytical models for predicting the Wsep were still generally limited to liquid-metal and oxide interfaces, and had to rely on some simple empirical correlations that described either the free energy of formation for oxide and liquid-metal, or the mixing enthalpies for respective oxide elements in metals [42]. Unfortunately, many of these models were not applicable to those systems comprised of non-oxidized ceramics interface. Furthermore, these models may only provide qualitative information for trends in adhesion and their generality. Even within the class of metal/oxide interfaces, it was still debatable because many of them only have those interfaces using the α-Al2O3 (Alumina) as oxides [43].

In recognition of limitations of analytic models, a number of ab-initio studies on metal/ceramic adhesion using density functional theory (DFT) have appeared during past two decades. These studies captured atomic-scale features at the interfacial areas and can reveal invaluable information using electronic structures. For example, one DFT research focused on the interfacial adhesion at several Al-Ceramic couples: Al/α-Al2O3, Al/WC, Al/VC, Al/VN, Al/CrN and Al/TiN [4447]. Among which, large alloy slabs obeying a unit stacking sequence along x- and y-directions, respectively, but up to 15 atomic layers along z-direction were adopted in constructing calculation geometries, so as to accurately evaluate the surface energies and avoid the non-convergence problem via the slab thickness. Superscripts on polar surface systems, α-Al2O3(0001)O, WC(0001)W, WC(0001)C, etc., indicated surface terminations. In addition, each non-stoichiometric slab must be constructed to allow for the identical termination of its both surfaces. Regarding their relevant geometric parameters, see [47]. In this research, although Eq. (10) defined the Wsep in terms of both surface and interface energies, these quantities do not necessarily play an equal role, that is, since interface energy γ is relatively small in those metal-ceramic systems, the strength of interfacial bonding relevant to the Wsep would depend to a large extent upon the reactivity of individual surfaces as reflected by their surface energies, see Table 1 [47].

If considering atomic relaxations into the calculation, at an Al/Alumina interface, both magnitude and rank ordering of adhesion energies for different stacking sequences are found to change to the unrelaxed results determined by a "universal binding energy relation (UBER, in this framework, both charge density and binding energy of surface via inter-particle separation were empirically described as exponential decaying functions)" [48], see Figure 1(a). Regarding the optimized interfacial stacking sequences, through analyzing electronic structures (charge density contours), it is found that, at the polar Al/α-Al2O3 and the polar Al/WC interfacial areas (more dangling bonds on isolated surfaces), metal atoms prefer to the sites which continue the natural stacking sequence in ceramic bulk across the interfacial area into the Al-slab host, see the right picture in Figure 1(b); while for non-polar Al/α-Al2O3Al, Al/VC, Al/VN, Al/CrN, Al/TiN and Al/WC systems, Al atoms prefer to the sits above metalloid atoms (O, C and N), see the left picture in Figure 1(b).


Figure 1.

(a) (Left) An optimized geometry for a fcc-Al/α-Al2O3 interface determined by the UBER calculations; (right) the final relaxed structure: Small spheres are Al atoms, and large spheres are O atoms. The view direction is along [12¯10]. (b) Side-view of electronic structures at an interface: Two slices through the electron localization function (ELF) for an fcc-Al/α-Al2O3 interface taken along (101¯0) and (112¯0) planes, indicating four of hexagonal close-packed oxygen layers in oxide (bottom) and all five atomic layers from an Al slab (top).

3.1.1. Effect of surface geometry to the Wsep

The Wsep calculations are listed in Table 2 [47], where readers may find out obvious differences in the Wsep values between polar and non-polar geometries; among which, three polar surfaces behave the largest Wsep values from 4.08 J⋅m–2 for the Al/WCW up to 9.73 J⋅m–2 for the Al/α-Al2O3O; for the remaining six non-polar interfaces, five have the Wsep values of 2.10 J⋅m–2 or less. Furthermore, values for three nitride-ceramics are among the smallest range and look relatively insensitive to the choice of metallic component, with the Wsep values falling within 1.45–1.73 J⋅m–2. Averagely, the Wsep values in polar systems are about 4.60 J⋅m–2 larger than in non-polar systems, consistent with their larger surface energies as shown in Table 1. It was also noted that, the predicted Wsep value of 1.06 J⋅m–2 for the Al/α-Al2O3Al system was in favorable agreement with the experimental value of 1.13 J⋅m–2 in [49].

SurfaceNo. of layersσ (J·m−2)

Table 1.

Calculated surface energies σi and slab-i dimensions (no. of layers) used for the Wsep calculations.

[i] - Superscripts give termination of those surfaces cleaved along a polar plane.

InterfaceOrientationPolarityStrain (%)γ (J·m−2)Wsep (J·m−2)

Table 2.

Interfacial orientation relationship, polarity (P = polar, NP = nonpolar), strain, interface energy (γ ), and Wsep.

Terminations of polar ceramic surfaces are indicated with a superscript. The Wsep values correspond to the optimized stacking sequences.

3.1.2. Effect of surface energies to the Wsep

To further verify whether the Wsep is mainly determined by surface energies or not, one Figure (see Figure 1 in [47]) plots the Wsep versus the [σAl+σCeramic]. A dashed-line with a unit slope is added to indicate data points lying in an idealized case: γ = 0. In this figure, the horizontal deviation of a given data point from the "γ = 0 – line" corresponds to the γ size. For three polar interfaces in the top region, data are bracketed with horizontal error bars, giving the range of possible σi listed in Table 1. The Wsep values in the whole range indicate a roughly linear trend with respect to the σi , all but one of data points fall to the line-right. One system falling to the line-left with a negative γ value and involving an α-Al2O3O surface, shows the largest surface energy and the Wsep value. This implies that the possibility of the fracture be within the Al-bulk rather than at the Al-Ceramic junction (interfacial area) due to a tensile stress. It also indicates that such an interface geometry is very firmly constructed.

With the exception of TiN, the rocksalt-structured carbides and nitrides (CrN, VN and VC) exhibit the largest dependence of the Wsep on their surface energies: These three systems have the smallest interface energies, only 0.03–0.18 J⋅m–2, and their Wsep values crowd in a nearly linear fashion (close to the γ = 0 – line) along the bottom left of the figure. Presumably, small γ values at these interfaces can be explained by similarities between bulk and interface bonding. However, since both of bulk components exhibit some extent to metallic bonding, for the Al/TiN system—a singularity with respect to other rocksalt ceramics, the γ value contributes larger, 0.62 J⋅m–2, which could be due to the relatively large interfacial strain (misfit of slab geometry) of 5.3%.

Correspondingly, Table 2 also illustrates that those Wsep values involving non-polar α-Al2O3 and WC (112¯0) systems are not solely depicted by their surface energies (i.e., the γ effects should be more important). For two α-Al2O3 interfaces, as with TiN, a possible explanation for this behavior could be the significant strain (4.9 and 5.3%). In addition, a main effect of ionic bonding in bulk α-Al2O3 differs substantially from that in Al-bulk. On the contrary for the WC system, the strain is very low (0.4%), and thus differences of bonding types in respective bulks are less appeared (WC is metallic). So it implements an incoherent, misfit geometry characterized by irregular bonding across the interface.

3.2. Conclusion

Comprehensively, the DFT calculations indicated that the surface energy played a dominant role in determining the adhesion energy (the Wsep) for systems described above, interface energy was generally much smaller in comparison. In most cases, particularly for those polar systems involving cubic oxides, carbides and nitrides (adhesives) at interfaces, the Wsep was described largely in terms of surface energy alone. Exceptions to this rule were mainly for systems with relatively large interfacial strain or incoherent interface bonding.

4. Ab-initio molecular dynamics study at an Al/Fe2O3 interface

In this section, without any loss of generality, ab-initio molecular dynamics (AIMD) based upon density functional theory (DFT-GGA) is employed to discuss the adhesive transfer of materials at an interface combining a clean Al(111) and a hematite [α-Fe2O3(0001)] slabs together, to obtain invaluable information on the dynamical adhesion. As a result, the compression of two slabs may form an amorphous alumina (Al2O3)-texture, which welds two slabs together at the interfacial area through a thermite reaction. Subsequently, a tensile process results in the removal of amorphous alumina-texture from the melting Al-slab. The overall reaction is highly exothermic to increase the temperature of the system.

4.1. Modeling geometry and simulation procedures

An AIMD study focusing on a clean Al(111) slab here involves a super-cell geometry with a periodic 4-layer Al(111)-slab (16 ions per layer) in the x- and the y-directions [50]. This orthorhombic super-cell has three orientated vectors: a[11¯0]=22aAl at the x-axis; b[112¯]=6aAl at the y-axis; and c[111]=40Å at the z-axis along with a vacuum distance of 24.0 Å in the z-direction, and one bottom layer fixed along the z-direction. Here aAl=4.053Å [51]. Similarly, simulations for a widely studied α-Fe2O3(0001) slab with oxygen-termination (more reactive than iron-termination to the Al-surface) are conducted using a periodic seven-layer rectangle slab with two oriented vectors: 5aFe2O3 at the x-axis and 152aFe2O3 at the y-axis, and a height of 8.40 Å vacuum space with free boundary along the z-direction (total of 36 Fe and 54 O ions). Here aFe2O3=5.031Å [52]. Comparing these two slabs, the Al-slab size is slightly larger than the Fe2O3-slab along x- and y-directions, respectively. So a super-cell fitting for the Al-slab is selected as the final frame to contain both of slabs in the simulation. In this super-cell, misfit percentages of each of two slab lengths along x- and y-directions are calculated at very small values: 1.9 and 1.9%, respectively (see Figure 2).

The simulation time step was selected as 0.001 ps (i.e., 1 ps = 10−12 s). The whole system was held at constant energy of AIMD simulation (no temperature control) when two slabs interacted. At the beginning of each simulation, top three layers of the α-Fe2O3-slab were constrained to move down along the z-direction at a constant velocity, VZ = – 4.00 Å/ps, corresponding to a displacement dZ = – 0.004 Å per time step. Since this speed, 400 m/s, was typical of real processing conditions, such shock loading may exceed the yield strength of the Al-slab, and the deformation of top layers in the Al-slab can be reasonably observed during contact simulation (Figure 2).


Figure 2.

Side-view of geometry at the Al(111)/α-Fe2O3(0001) interface.

Before carrying out the compressive simulation, a favorable interface geometry between two slabs is optimized using the energy minimization in the DFT. The result indicates that the final spacing between two slabs would be about 1.88 Å to form Al–O bonds at the interfacial area, which is in very agreement with those observed in experiments [53] (see Figure 3). Table 3 also lists some relevant bond lengths [50].


Figure 3.

The final Al(111)/α-Fe2O3(0001) interfacial geometry optimized by the DFT.

Bond pairAl–OFe–O
Bond length (Å)CalculationExperiment [a]CalculationExperiment [b]

Table 3.

Relevant bond lengths at the interfacial area [unit: Å].

[a] Ref. [55].

[b] Ref. [56].

Therefore in the following AIMD simulation, the initial spacing between two slabs was selected as 2.50 Å, which was far enough to avoid the significant forming of chemical bonds in the R–L region (the interfacial area) when the simulation started. As two slabs moved closer, chemical bonds would begin forming in the spacing, and the interfacial reactions may occur. In Figure 4(a), since top three layers (outside the R–L region) in the α-Fe2O3(0001)-slab may move down at a constant velocity (VZ = − 4.0 Å/ps) as a rigid block, its final downward displacement is chosen as 0.62 Å, so as to keep the final spacing between two slabs around an optimized value: 1.88 Å.


Figure 4.

Main sequence of compressing the Al(111)/α-Fe2O3(0001) interface without lubricants.

4.2. Results and discussion

4.2.1. Compressive process

Figure 4 shows the compressing sequence of the Al(111) and the α-Fe2O3(0001) slabs. A fixed velocity, VZ = − 4.00Å/ps, is applied onto top three layers in the α-Fe2O3(0001)-slab along the z-direction during the simulation [see Figure 4(a)].

In this figure, during the first 80 time steps, the α-Fe2O3(0001)-slab is moving downward the Al(111)-slab without forming Al–O bonds, see Figure 4(a); around the 80th time step, some Al–O bonds begin forming (with Al–O bond length about 1.86~1.93 Å) [see Figure 4(b)]; around the 140th time step, most of Al–O bonds in the R–L region are formed [see Figure 4(c)]. Local temperature in this region increases to about 870 K. During this evolution, an amorphous Al2O3-texture is gradually formed due to the further reaction of oxygen-ions with the Al-surface; around the 390th time step, this reaction continues at the uppermost point [see Figure 4(d)]. Local temperature in the R–L region continues rising to about 1130 K; by the 620th time step, some new layers made up of amorphous alumina-texture are completely constructed in the R–L region [see Figure 4(e)]. Local temperature in this region increases to about 1200 K.


Figure 5.

The stress perpendicular to the interfacial area as two slabs compress together.

Figure 5 indicates the distribution of the stress normal to the interfacial area as two slabs compress together. Initially, the stress becomes negative (attractive between two slabs) as Al–O bonds begin forming in the R–L region (points a–c); when two slabs are moving closer, the stress tends to the positive values (repulsive between two slabs) because most of Al and O atoms at the interfacial area are forced to come closely, leading to a nearly formation of amorphous alumina-texture with thermal (exothermic) expansion; after two slabs totally react, the stress falls down to a near zero point (point e) from the highest positive one (point d) as an optimal bond spacing between two slabs is reached; and then, the stress increases again (after point e) when two slabs continue compressing. This figure also represents a polynomial (dashed line) fit to the stress curve, which may clearly regress the real stress trend through its oscillation. Please note that the stress trend just depicts the starting stage of a whole P-h curve as stated in [54]. However, since two slabs here are too thin and too small in dimension, and at too short time steps throughout the whole AIMD simulation, it is hardly to rough out the whole P-h curve fulfilled by classical molecular dynamics on very large size models at very long time steps. Even if, this stress curve may still provide some invaluable information at the initial stage of nanoindentation process.

Figure 6 shows the number change of Al–O and Fe–O bonds versus normal distance at the interfacial area in Figure 4(a) during the compressive process. In this figure, Al–O bonds in the R–L region increase substantially while Fe–O bonds in this region decrease gradually. The reason is that, when two slabs compress together, the dangling bonds on Al and O atoms at the interfacial area would bond together. As the compression continues, more Al–O bonds are forming as some Fe–O bonds are breaking. Nevertheless, only a very porous interface is formed during two slabs compression, which is far from the ideally compressed status. Therefore, the thermite reaction at the interfacial area just partially occurs, only a few Fe–O bonds are breaking in this region.


Figure 6.

The number change of Al–O and Fe–O bonds in the R–L region during slab compression.

4.2.2. Tensile process

After the compression, an effect of a tensile strain acting on the compressed α-Fe2O3(0001)/Al(111) slabs is immediately conducted: Top three layers of the α-Fe2O3(0001)-slab in Figure 4(e) is moving up as a rigid block at a constant velocity, VZ = + 4.00 Å/ps, reversing the downward direction of previous compression.

Figure 7 shows the stretching sequence of this process: (i) when top three layers of the α-Fe2O3(0001)-slab move up about 1.30 Å from the starting point in Figure 4(e), the real tensile process at the interfacial area begins (part a). (ii) Then next, as two slabs separate, further deformation at the interfacial area and thermite reaction will continue, more atomic mixing in this area oxidizes the Al-surface (part b). (iii) During transitional stages (parts c and d), the whole interfacial area begins breaking apart, mainly due to the breaking of some more Fe–O and Al–Al bonds. Along with a small increase of Al–O bonds, amorphous alumina-texture is continuously formed, and an adhesive metal transfer begins. (iv) Finally, surface on the Al-slab is totally broken: an adhesive metal transfer has completely occurred, amorphous alumina-texture is removed from the remnant surface of the Al-slab, and bonds onto the bottom of the α-Fe2O3(0001)-slab, indicating the final apart step of adhesive transfer of materials. In addition, a fresh disordered Al-surface is exposed on the residual Al-slab after the tensile process finishes (part e).


Figure 7.

The main sequence of stretching the Al(111)/α-Fe2O3(0001) slabs.


Figure 8.

The number change of Al–O, Fe–O and Al–Al bonds in the new R–L region during the slab tension.

Figure 8 shows the number change of Al–O, Fe–O and Al–Al bonds in a new larger R–L region as shown in Figure 7(a) during the tensile process. In this figure, at initial stage, a few more Al–O bonds are forming and a few more Fe–O bonds are breaking, which represents the continuous thermite reaction in the R-L region as the Fe2O3-slab is being pulled upward. Also, there is a breaking of approximately 10–12 Al–Al bonds in top two layers of the Al-slab. When the Fe2O3-slab is pulled upward about 8.0 Å, there is a slight increase in Al–O bonds but a slight decrease in Fe–O bonds. In this case, the interfacial area is totally separated. By the end of the tensile process, the Fe–O coordination is approximately 1:1, while the Al–O one is approximately 2:3. Thus, this reaction could be approximately 2Al + 3Fe2O3 → 6FeO + Al2O3, instead of a more favorable thermite reaction: 2Al + 3FeO → 2Al2O3 + 3Fe, which is mainly due to very limited reaction time. Local temperature in the residual Al-slab is about 810 K, while that in the removed slab (Fe2O3/FeO+Al2O3) is about 3030 K, implying that such a soft slab may bond onto a hard-roller surface as a protective thin film resistant to the elevated temperature.

Comprehensively, the AIMD simulations reveal that the tensile stress at the Al(111)/α-Fe2O3(0001) interface may remove the Al2O3 fragments from the Al-slab surface after a steel roller (covered with an iron oxide) contacts the slab under an external strain load, which could indicate a typical process of the adhesive transfer of materials at interfaces.

4.3. Conclusion

The AIMD simulations of adhesive transfer of materials at an interfacial area without lubricants were performed for compressive and tensile processes between the Al(111) and the α-Fe2O3(0001) slabs. The conclusions can be drawn: (1) Contact of the Fe2O3-slab with the Al-slab may result in a strong exothermic reaction to form an amorphous alumina-texture, which welded two slabs together. (2) Separation of two slabs may have the local temperature at the interfacial area increased significantly, leading to the melting of the Al-surface and easy removal of amorphous Al2O3-texture from the Al-slab surface. (3) The initial stage of the P-h curve can be drawn throughout the AIMD simulation to provide the invaluable information, which was hardly observed throughout classical molecular dynamics simulation. (4) During the AIMD simulation, the whole reaction process at the interfacial area was 2Al + 3Fe2O3 → 6FeO + Al2O3, that is, due to the very short simulation time, such a thermite reaction was incomplete.

5. Summary

In this chapter, atomistic simulations of adhesive behaviors at several Al/Ceramic/Iron oxide interfaces are comprehensively discussed using density functional theory and ab-initio molecular dynamics, to illuminate a reasonable way for people’s further understanding to this fascinating field. Some main points in the discussion include, without any loss of generality:

  1. Main phenomenon in the Al/Ceramic adhesion indicated: In most cases, surface energy would play a dominant role in determining adhesion energy (Wsep) as interface energy was generally much smaller in comparison if ignoring effects of large stain and incoherent bonding at the interfacial area between two contact slabs, particularly for those polar systems involving cubic oxides, carbides and nitrides to form very solid thin films (adhesives) at interfaces. Apparently, the more the Wsep depended upon the surface energy, the more steady and tenacious the adhesives (thin films) bind the interfaces. Therefore, knowledge of surface energies alone could serve as a reasonable starting point for estimating interfacial strength.

  2. Common phenomena in the adhesive transfer of materials (phase transformation of indented materials) included: (a) Contact of the Fe2O3-slab with the Al-slab may result in a strong exothermic reaction to form an amorphous alumina-texture, which welded two slabs together. (b) Separation of two slabs led to the melting of the Al-surface and easy removal of amorphous Al2O3-texture from the Al-slab surface. This texture may bond onto a hard-roller surface as a protective thin film (soft adhesives) resistant to the elevated temperature.


This work was supported by the National Science Foundation (Grant no.: DMR9619353), USA; the Hebei Province Science and Technology Support Program (Grant no.: 15961006D), the Start-up Fund of Doctoral Research at North China Institute of Aerospace Engineering (Grant no.: BKY201405), and the Opening Fund of State Key Laboratory of Nonlinear Mechanics at Institute of Mechanics in Chinese Academy of Sciences (Grant no.: LNM201503), China.


1 - Doerner MF, Nix WD. A method for interpreting the data from depth-sensing indentation instruments. Journal of Materials Research. 1986; 1: 601-609. DOI:
2 - Oliver WC, Pharr GM. An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments. Journal of Materials Research. 1992; 7: 1564-1583. DOI:
3 - Gerberich WW, Nelson JC, Lilleodden1 ET, Anderson E, Wyrobek JT. Indentation induced dislocation nucleation: The initial yield point. Acta Materialia. 1992; 44: 3585-3598. DOI: 10.1016/1359-6454(96)00010-9
4 - Gouldstone A, Koh HJ, Zeng KY. Discrete and continuous deformation during nanoindentation of thin films. Acta Materialia. 2000; 48: 2277-2295. DOI: 10.1016/S1359-6454(00)00009-4
5 - Page TF, Oliver WC, McHargue CJ. The deformation behavior of ceramic crystals subjected to very low load (nano) indentations. Journal of Materials Research. 1992; 7: 450-473. DOI:
6 - Pollock HM, Blau PJ. Friction, Lubrication, and Wear Technology-ASM Metals Handbook. Vol.18. ASM International; 1992. 419p. DOI: 10.1007/s11191-012-9502-4
7 - Shenoy VB, Miller R, Tadmor EB, Phillips R, Ortiz M. Quasi-continuum models of interfacial structure and deformation. Physical Review Letters. 1998; 80: 742-745. DOI:
8 - Ortiz M, Cuitino AM, Knap J, Koslowski M. Mixed atomistic-continuum models of material behavior: the art of transcending atomistics and informing continua. MRS Bulletin. 2001; 26: 216-221. DOI:
9 - Born M, Oppenheimer R. Zur Quantentheorie der Molekeln (On the Quantum Theory of Molecules). Annalen der Physik (in Germany). 1927; 389: 457-484. DOI: 10.1002/andp.19273892002
10 - Hohenberg P, Kohn W. Inhomogeneous electron gas. Physical Review. 1964; 136: B864-871. DOI: 10.1103/PhysRev.136.B864
11 - Kohn W, Sham LJ. Self-consistent equations including exchange and correlation effects. Physical Review. 1965; 140: A1133-1138. DOI: 10.1103/PhysRev.140.A1133
12 - Kresse G, Joubert D. From ultrasoft pseudopotentials to the projector augmented-wave method. Physical Review B. 1999; 59: 1758-1775. DOI:
13 - Jones RO, Gunnarsson O. The density functional formalism, its applications and prospects. Reviews of Modern Physics. 1989; 61: 689-746. DOI:
14 - Perdew, JP. Accurate density functional for the energy: real-space cutoff of the gradient expansion for the exchange hole. Physical Review Letters. 1985; 55: 1665-1668. DOI:
15 - Perdew JP, Chevary JA, Fiolhais C. Atoms, molecules, solids, and surfaces: applications of the generalized gradient approximation for exchange and correlation. Physical Review B. 1992; 46: 6671-6687. DOI:
16 - Perdew, JP. Generalized gradient approximations for exchange and correlation: a look backward and forward. Physica B. 1991; 172: 1-6. DOI: 10.1016/0921-4526(91)90409-8
17 - Perdew JP, Burke K, Ernzerhof M. ERRATA: generalized gradient approximation made simple. Physical Review Letters. 1996; 77: 3865-3868. DOI:]. Physical Review Letters. 1997; 78: 1396, DOI:
18 - Chadi DL, Cohen ML. Special points in the Brillouin-zone. Physical Review B. 1973; 8: 5747-5753. DOI:
19 - Monkhorst HJ, Pack JD. Special points for Brillouin-zone integrations. Physical Review B. 1976; 13: 5188-5192. DOI:
20 - Bachelet GB, Hamann DR, Schlüter M. Pseudopotentials that work: from H to Pu. Physical Review B. 1982; 26: 4199-4228. DOI:
21 - Vanderbilt D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Physical Review B. 1990; 41: 7892(R)-7895. DOI:
22 - Blöchl PE. Projector augmented-wave method. Physical Review B. 1994; 50: 17953-17979. DOI:
23 - Kresse G, Hafner J. Norm-conserving and ultrasoft pseudopotentials for first-row and transition elements. Journal of Physics: Condensed Matter. 1994; 6: 8245-8258. DOI: 10.1088/0953-8984/6/40/015
24 - Alder BJ, Wainwright TE. Phase transition for a hard sphere system. Journal of Chemical Physics. 1957; 27: 1208. DOI:; ibid. Studies in Molecular Dynamics. I. General Method. Journal of Chemical Physics. 1959; 31: 459-466, DOI:
25 - Verlet L. Computer "experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Physical Review. 1967; 159: 98-103, DOI:
26 - Beeman D. Some multistep methods for use in molecular dynamics calculations. Journal of Computational Physics. 1976; 20: 130-139. DOI: 10.1016/0021-9991(76)90059-0
27 - Gear CW. Numerical Initial Value Problems in Ordinary Differential Equations. 1st ed. Prentice-Hall, Upper Saddle River; 1971. 253p. ISBN-13: 978-0136266068
28 - Heermann DW. Computer Simulation Method in Theoretical Physics. 2nd ed. Springer Verlag, Berlin; 1990. 148p. ISBN-13: 978-3540522102
29 - Nosé S. A molecular dynamics method for simulations in the canonical ensemble. Molecular Physics. 1984; 52: 255-268. DOI: 10.1080/00268978400101201
30 - Anderson HC. Molecular dynamics simulations at constant pressure and/or temperature. Journal of Chemical Physics. 1980; 72: 2384-2393. DOI:
31 - Parrinello M, Rahman A. Polymorphic transitions in single crystals: a new molecular dynamics method. Journal of Applied Physics. 1981; 52: 7182-7190. DOI:
32 - Hoover, WG. Computational Statistical Mechanics. 1st ed. Elsevier Science, Oxford; 1991. 330p. ISBN-13: 978-0444564238
33 - Kaplan IG. Intermolecular Interactions-Physical Picture, Computational methods and Model Potentials. 1st ed. Wiley; 2006. 380p. ISBN-13: 978-0470863329
34 - Payne MC, Allan DC, Arias TA, Joannopoulos JD. Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients. Reviews of Modern Physics. 1992; 64: 1045-1098. DOI:
35 - Marx D, Hütter J. Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods. 1st ed. Cambridge; 2009. 577p. ISBN-13: 978-1107663534
36 - Car R, Parrinello M. Unified approach for molecular dynamics and density-functional theory. Physical Review Letters. 1985; 55: 2471-2474. DOI:
37 - VASP the GUIDE [Internet]. 2016. Available from:
38 - Henkelman G, Jónsson H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. Journal of Chemical Physics. 2000; 113: 9978-9985. DOI:
39 - Finnis, MW. The theory of metal-ceramic interfaces. Journal of Physics: Condensed Matter. 1996; 8: 5811-5836. DOI: 10.1088/0953-8984/8/32/003
40 - Humenik M, Kingery WD. Metal-ceramic interactions: III, surface tension and wettability of metal-ceramic systems. Journal of the American Ceramic Society. 1954, 37: 18-23. DOI: DOI: 10.1111/j.1151-2916.1954.tb13972.x
41 - Howe JM. Bonding, structure, and properties of metal/ceramic interfaces: Part 1 Chemical bonding, chemical reaction, and interfacial structure. International Materials Reviews. 1993; 38: 233-256. DOI:
42 - Eustathopoulos N, Chatain D, Coudurier L. Wetting and interfacial chemistry in liquid metal-ceramic systems. Materials Science and Engineering A. 1991; 135: 83-88, DOI: 10.1016/0921-5093(91)90541-T
43 - Barrera RG, Duke CB. Dielectric continuum theory of the electronic structure of interfaces. Physical Review B. 1976; 13: 4477-4489. DOI:
44 - Siegel DJ, Hector LG Jr, Adams JB. Adhesion, atomic structure, and bonding at the Al(111)/α−Al2O3(0001) interface: a first principles study. Physical Review B. 2002; 65: 085415. DOI:
45 - Siegel DJ, Hector LG Jr, Adams JB. Adhesion, stability, and bonding at metal/metal-carbide interfaces: Al/WC. Surface Science. 2002; 498: 321-336. DOI: 10.1016/S0039-6028(01)01811-8
46 - Siegel DJ, Hector LG Jr, Adams JB. First-principles study of metal–carbide/nitride adhesion: Al/VC vs. Al/VN. Acta Materialia. 2002; 50: 619-631. DOI: 10.1016/S1359-6454(01)00361-5
47 - Siegel DJ, Hector LG Jr, Adams JB. Ab initio study of Al-ceramic interfacial adhesion. Physical Review B. 2003; 67: 092105. DOI:
48 - Banerjea A, Smith R. Origins of the universal binding-energy relation. Physical Review B. 1988; 37: 6632-6645. DOI:
49 - Lipkin DM, Israelachvili JN, Clarke DR. Estimating the metal-ceramic van der Waals adhesion energy. Philosophical Magazine A. 1997; 76: 715-728. DOI: 10.1080/01418619708214205
50 - Zhong J, Adams JB. Adhesive metal transfer at the Al(1 1 1)/α-Fe2O3(0001) interface: a study with ab initio molecular dynamics. Modelling and Simulation in Materials Science and Engineering. 2008; 16: 085001. DOI: 10.1088/0965-0393/16/8/085001.
51 - Zhong J, Adams JB. Adsorption and decomposition pathways of vinyl phosphonic and ethanoic acids on the Al(111) surface:  a density functional analysis. Journal of Physical Chemistry C. 2007; 111: 7366-7375. DOI: 10.1021/jp0667487
52 - Yamamoto N. The shift of the spin flip temperature of α-Fe2O3 fine particles. Journal of the Physical Society of Japan. 1968; 24: 23-28. DOI:
53 - Coast R, Pikus M, Henriksen NP, Nitowski AG. A vibrational spectroscopic comparison of vinyl-triethoxysilane and vinyl-phosphonic acid adsorbed on oxidized aluminum. Journal of Adhesion Science and Technology. 1996; 10: 101-121. DOI: 10.1163/156856196X00805
54 - Szlufarska I. Atomistic simulations of nanoindentation. Materials Today. 2006; 9: 42-50. DOI: 10.1016/S1369-7021(06)71496-1
55 - Eng PJ, Trainor TP, Brown Jr GE, Waychunas GA, Newville M, Sutton SR, Rivers ML. Structure of the Hydrated α-Al2O3(0001) Surface. Science. 2000; 288: 1029–1033. DOI: 10.1126/science.288.5468.1029.
56 - Hay MT, Geib SJ. Tetra-butyl-ammonium isobutyl silsesquioxane monochloro-ferrate(III). Acta Crystallographica Section E. 2007; E63: m445-m446. DOI: 10.1107/S1600536807000360.