Effects of van der Waals Interactions in the Adsorption of Isooctane and Ethanol on Fe(100) Surfaces.

van der Waals (vdW) forces play a fundamental role in the structure and behavior of diverse systems. Because of development of functionals that include nonlocal correlation, it is possible to study the effects of vdW interactions in systems of industrial and tribological interest. Here we simulated within the framework of density functional theory (DFT) the adsorption of isooctane (2,2,4-trimethylpentane) and ethanol on an Fe(100) surface, employing various exchange-correlation functionals to take vdW forces into account. In particular, this paper discusses the effect of vdW forces on the magnitude of adsorption energies, equilibrium geometries, and their role in the binding mechanism. According to our calculations, vdW interactions increase the adsorption energies and reduce the equilibrium distances. Nevertheless, they do not influence the spatial configuration of the adsorbed molecules. Their effect on the electronic density is a nonisotropic, delocalized accumulation of charge between the molecule and the slab. In conclusion, vdW forces are essential for the adsorption of isooctane and ethanol on a bcc Fe(100) surface.

Introduction

With a universal presence in all molecules and solids, van der Waals forces (vdW) are important interactions to consider in the study of matter.[1−3] Although weak in comparison to chemical bonds, their long-range nature and their collective effect play a decisive role in the structure of molecules and their interaction with surfaces.[4−7] Nevertheless, their description from first principles has proven to be challenging. Most of the common theoretical methods, such as the density functional theory (DFT), fail to describe them properly. Typically, the commonly used generalized gradient approximation (GGA) underestimates the binding energies and overestimates the equilibrium distances in various physisorbed systems, while the ones calculated within the local density approximation (LDA) are closer to the experiment. However, this better agreement is fortuitous, since the exponential decline of the LDA interaction cannot account for the polynomial long-range behavior of the vdW interactions.[8,9]

Over the past decade, several approximations have been proposed[10−39] to take dispersion into account. The most sophisticated methods aim to treat dispersion beyond the pairwise approximation by considering collective excitations. The many-body dispersion approach[24,25] is one of these models that uses coupled dipoles. The dispersion interaction is obtained from shifts in the frequencies of harmonic oscillators that occupy the atomic positions, as the interaction between them is activated. Although promising, it is challenging to get accurate relations between atoms and oscillator models.[40] Another model that includes the vdW energy accurately is the random phase approximation,[27,29−31] combined with the adiabatic connection and fluctuation dissipation theorem.[26,28] The computational cost of this method limits it at present, however, to small systems commonly used as benchmarks. Other models have been proposed that consider dispersion to be pairwise additive, but they are otherwise independent of any external input parameters, such as VV10,[38] the local response dispersion approach,[13,14] and the van der Waals density functional (vdW-DF).[10] Klimeš and co-workers[41,42] proposed a series of optimized functionals within the framework of DFT based on the vdW-DF, which have proven to be among the most accurate in this class of models. Furthermore, the vdW-DF and its optimized versions have been implemented in extensively distributed DFT codes and applied to a wide range of systems including dimers,[41] soft layered materials,[9] organic molecules adsorbed on graphite,[43,44] as well as graphene and noble gases on metals,[8,45,46] among others. These works showed the advantages and disadvantages of these methods and consequently increased their reliability on subsequent applications.

Of special interest for technological purposes are the interactions between organic molecules and metallic surfaces. The effect of nonlocal forces has been previously investigated in the adsorption process occurring in promising candidates for optoelectronic devices, including thiophene on Cu(110)[47] and benzene, along with related compounds, on several transition metals.[48,49] The influence of dispersion forces has also been studied in the adsorption of n-butane on copper and gold surfaces,[50] which are considered prototypes of typical weak physisorption systems. However, in many industrial applications the metallic surfaces are iron-based alloys. Moreover, the involved molecules may have larger size, be branched, or contain other types of functional groups, as in the case of fuels and lubricants as well as their base fluids and additives. The usefulness of these compounds strongly depends on the interacting forces at the molecule–surface interface, since they play an important role in the tribological behavior; i.e., low friction and low wear are the most desired requirements for their application.

Recent work in the field of fuel tribology using molecular dynamics (MD) simulations covered the stability of monolayers of stearic acid adsorbed onto nanorough iron surfaces during shearing[51,52] or the frictional performance of fuel additives.[53,54] On the microscale, a joint numerical approach combining the finite-element method and the boundary-element method (FEM-BEM) has been employed to simulate wear processes in a diesel fuel lubricated sliding contact.[55] However, none of these nano- and microtribological simulations using MD or FEM-BEM explicitly considered the adsorption process as an initial step toward the lubricity of surface-active species. Moreover, the knowledge of the adsorption behavior of ethanol is a crucial preliminary step for a better description of the effectiveness of additives in fuels containing biocomponents, i.e., substances with considerably higher polarity in comparison with conventional aliphatic and aromatic fuel components.

Last year, Tereshchuk and Da Silva[56] studied the adsorption of ethanol and water on several transition-metal surfaces, including Fe(110). Although the reported adsorption energies calculated with vdW corrections are notably different from the corresponding GGA results, the vdW correction applied in that study belongs to a class of empirical approaches that are computationally efficient but overall less accurate than other available approaches to take vdW into account.[40] Thus, a comparison with a method independent of external input parameters is desirable to validate the accuracy of the calculations. Furthermore, a detailed insight into the adsorption mechanism and the influence of the surface termination is required for a better description of these systems. To study the effects of vdW interactions in the adsorption process of molecules and surfaces of industrial and tribological interest, we simulated via first-principles calculations the adsorption of isooctane (2,2,4-trimethylpentane) and ethanol on a bcc Fe(100) surface employing various GGA and vdW functionals. These adsorbates were selected because isooctane is a representative aliphatic gasoline compound and ethanol can be found (besides gasoline) in many relevant products. Additionally, some properties differ remarkably in these molecules, such as the chemical polarity, polarizability, and chemical reactivity, and it is therefore possible to investigate their influence on the adsorption process. This paper describes the magnitude of adsorption energies calculated with a rigorous treatment of nonlocal interactions as well as the influence of vdW forces on equilibrium geometries and their role in the binding mechanism.

Computational Details

To study the adsorption energies and equilibrium distances, we carried out spin-polarized first-principles calculations within the framework of DFT.[57] The Vienna ab initio simulation package (VASP)[58−63] was used to perform the required computations. VASP produces an iterative solution of the Kohn–Sham equations within a plane-wave basis, employing periodic boundary conditions. The projector augmented wave (PAW) method[64] was applied to describe the interaction between the core and the valence electrons.

We considered various approximations to the exchange and correlation functionals. Initially, we tested two functionals constructed in the generalized gradient approximation (GGA): the Perdew, Burke, and Ernzerhof (PBE)[65,66] and its revised version (revPBE).[67] These functionals differ only in one parameter of the exchange term, κ, changed from 0.804 in PBE to 1.245 in revPBE, and both have been extensively applied in physics and chemistry. In a second step, we employed the vdW density functional (vdW-DF) by Dion et al.,[10]where ErevPBE(x) [n] is the exchange energy obtained with the revPBE functional, ELDA(c) [n] is an LDA correlation and Enl(c) [n] is a nonlocal correlation term that approximates the vdW interactions. We also considered the optimized Becke86[68] van der Waals (optB86b-vdW) functional introduced by Klimeš et al.,[41,42]where a reparametrized version of the Becke86 exchange functional replaces the revPBE exchange used in eq 1. Among the two functionals described in eqs 1 and 2, the optB86b-vdW is generally more accurate,[42] and the results obtained with it should be preferred as reference. However, we apply the vdW-DF to analyze the effects of nonlocal correlation, since we can subtract the contributions of the other terms by introducing a revPBE+LDA functional,which involves terms that have been tested and applied in a wide range of materials and whose behavior, in contrast to the exchange term in eq 2, is well-known.

To model our systems, we constructed a supercell consisting of a body-centered cubic (bcc) iron slab and a molecule placed on top of it. Ten layers of atoms were included in the slab, each one containing 25 iron atoms. The vacuum spacing in the z-direction of repeated cells was 29.73 Å. The z axis is parallel to the long axis of the simulation cell and starts at the bottom layer of the iron slab. This setup accurately models a bcc Fe(100) surface and avoids molecule–molecule and slab–slab interactions, as well as the need to include dipole corrections. The lattice parameter of 2.83 Å used to construct the iron slab was obtained from a bulk calculation for bcc iron using the PBE functional, where the calculation parameters were chosen to keep the accuracy consistent with the rest of the computations. This value is in good agreement with the experimental lattice constant of 2.86 Å.[69]

To ensure sufficiently accurate total energies and forces, we carefully selected and tested our calculation parameters. A tight convergence criterion of 10–6 eV on the total energy in the self-consistency cycle was used, and a cutoff energy of 400 eV was applied for the plane-wave basis set. The k-space integrations were performed using a 2 × 2 × 1 Monkhorst–Pack mesh,[70,71] whereas the tetrahedron method with Blöchl corrections was employed for the static calculations and a Gaussian smearing with a width of 0.2 eV for the relaxations. The conjugate gradient algorithm was used to relax the structures, allowing the ions to move until an energy convergence criterion of 10–5 eV was fulfilled.

To find the equilibrium structures, we calculated the total energies of varied spatial configurations followed by a relaxation of the most stable one. Various adsorption sites, slab–molecule distances, and molecule orientations were analyzed to obtain a better starting guess for the relaxations. During relaxations, the ions in the top four layers of the slab and the ones constituting the molecule were allowed to move in all directions, while the atoms in the remaining layers of the slab were kept fixed to simulate the bulk properties.

We calculated the potential energy curves by varying the distance between the molecule and iron slab and subsequently computing the total energy of the resulting system via static calculations. The starting geometry for the displacements was the equilibrium structure relaxed with the corresponding functional. The separation between the molecule and the slab (d) was defined as the vertical distance between the atom with the lowest z coordinate in the molecule and the iron atom closest in distance to it. The equilibrium distance was considered to be the separation (as defined above) between the molecule and the slab that minimizes the energy of the system. We stress that this distance is obtained from a finite set of energy points calculated for various values of d, and therefore, its accuracy depends on the step size between those values.

Results and Discussion

Neither the orientation nor the adsorption site of isooctane and ethanol is influenced by the nonlocal correlation. Initially, the preferred adsorption configurations were calculated within the GGA aproximation. After inclusion of nonlocal interactions in our calculations via the optB86b-vdW functional, a second relaxation did not significantly change the geometry of this configuration. The largest variations in bond lengths were of the order of 10–3 Å, and in angles, the variations were of the order of 10–3 degree. In the most energetically favorable orientation, the carbon atoms of the isooctane molecule are close to the top and hollow sites of the iron slab (Figure 1). The energies of several other orientations differ by only around 20 meV per supercell, and for this reason no orientation is particularly favored at room temperature, where thermal energy kBT = 25 meV. Similarly, in the adsorption of ethanol several orientations of the molecule are possible. In this case, however, the hydroxyl group always orients itself toward the slab in all low-energy configurations (Figure 2). The plotted geometry differs from the one reported by Tereshchuk and Da Silva[56] for ethanol adsorbed on the Fe(110) surface only by the C–C bond, which is almost perpendicular to the surface. Our calculated small energy difference of 3 meV between these two configurations indicates that both can coexist at room temperature and that the surface termination does not influence the orientation of the adsorbed molecule. The energy hierarchy of all PBE structures did not change when it was recalculated with optB86b-vdW to take nonlocal forces into account.

Figure 1

Equilibrium adsorption geometry of an isooctane molecule on a bcc Fe(100) surface. The top, bridge, and hollow positions are indicated by a circle, triangle, and square, respectively.

Figure 2

Equilibrium adsorption geometry of an ethanol molecule on a bcc Fe(100) surface. The top, bridge, and hollow positions are indicated by a circle, triangle, and square, respectively.

Nonlocal interactions increase the adsorption energy of isooctane on a bcc Fe(100) surface and reduce the equilibrium distance (Table 1). The absolute value of the adsorption energy calculated with the optB86b-vdW functional is more than 12 times larger than the one calculated with the PBE functional, while the distance between the slab and the molecule, when in equilibrium, is reduced by 1.00 Å (Figure 3). To investigate the contribution of the nonlocal correlation to this increment, we compared the adsorption energies calculated with the vdW-DF and the revPBE+LDA functionals. The resulting energy difference is then 418 meV (Figure 3). Since this difference accounts for 95.8% of the adsorption energy calculated with the vdW-DF functional, the variation can now be unambiguously attributed to the dispersion forces. By the same reasoning, the decrease in the binding distance of 1.50 Å can also be associated with nonlocal interactions.

Table 1

Adsorption Energies and Equilibrium Distances Calculated with Various Exchange–Correlation Functionals

	ads energy (meV)		eq distance (Å)a
functional	isooctane	ethanol	isooctane	ethanol
PBE	47	384	3.00	2.00
revPBE+LDA	18	26	4.00	4.00
revPBE		144		2.25
optB86b	585	795	2.00	2.00
vdW-DF	436	517	2.50	2.25

As previously mentioned, the equilibrium distance is determined from a finite set of structures with fixed values of d, with a step size of 0.25 Å or smaller.

Figure 3

Calculated adsorption energy of an isooctane molecule on the bcc Fe(100) surface as a function of the vertical distance for the PBE, revPBE+LDA, vdW-DF, and optB86b-vdW exchange–correlation potentials.

As a result of the dispersion forces, the nonlocal correlation between electrons induces the change in the adsorption energy. The equilibrium distance is mainly determined by a balance between the long-range attractive vdW forces and the short-range Pauli repulsion. When the isooctane molecule approaches the iron slab, the Pauli repulsion causes a redistribution of the charge density, particularly among the d-orbitals of the iron ions. No charge is transferred between the molecule and the iron slab during this process. The overlapping between the wave functions of the molecule and the slab accounts for the movement of the electrons to higher energy states, increasing the total energy of the system. This effect is weaker when a proper description of nonlocal interactions is considered because the nonlocal correlations reduce the electron–electron repulsion (Figure 4). This allows the isooctane molecule to reach a shorter equilibrium distance, where the magnitude of the attractive forces is larger and, consequently, the binding energy increases. The calculated equilibrium distance is also affected by the choice of the vdW density functional, since the Pauli repulsion gives rise to exchange interactions and these functionals differ in the description of the exchange energy. For instance, the difference between the binding distance calculated with the optB86b-vdW functional and the one calculated with the vdW-DF is 0.50 Å.

Figure 4

Charge density difference (ρdiff) of isooctane adsorbed on the bcc Fe(100) surface at d = 2.00 Å. The charge density difference is defined as ρdiff = ρ – (ρisooctane + ρFe(100)) where ρ denotes the charge density of isooctane adsorbed on Fe(100), while ρisooctane and ρFe(100) represent the charge densities of the isolated molecule and the clean Fe(100) surface, respectively. The charge density difference is plotted in a plane perpendicular to the surface for values between −5 × 10–4 (solid blue, deficit) and 5 × 10–4 (solid red, accumulation) electrons/Å3.

As in the case of isooctane, the nonlocal correlation enhances the binding energy between the adsorbed ethanol molecule and the bcc Fe(100) surface (Table 1). The adsorption energy calculated with the optB86b-vdW functional is only 2 times larger than the one calculated with the PBE functional, and the equilibrium separation is 2.00 Å in both cases (Figure 8). This difference is remarkably smaller than the one in the adsorption of isooctane and can be caused either by the involvement of other forces in the binding mechanism or by an artifact, as in the LDA functional. To investigate these possibilities, the effect of the long-range interactions is extracted by comparing calculations with the vdW-DF functional to calculations with revPBE+LDA, as previously described. The contribution to the adsorption energy which can be attributed to the dispersion forces is 95% (Figure 8). Although this shows a significant contribution of the nonlocal interactions to the binding mechanism, the difference of 278 meV between the adsorption energies calculated with the vdW-DF and the optB86b-vdW functional needs to be investigated before discarding the contribution of other forces to the adsorption energy.

Figure 8

Calculated adsorption energy of an ethanol molecule on the bcc Fe(100) surface as a function of the vertical distance for the PBE, revPBE, revPBE+LDA, vdW-DF, and optB86b-vdW exchange–correlation potentials.

The nonlocal interactions contribute to the binding process of isooctane by inducing an accumulation of charge between the iron surface and the adsorbed molecule. Thonhauser et al.[72] showed the nature of the vdW bonds in the argon dimer by comparing the induced electron density between two calculations that differ in the inclusion of a term corresponding to the nonlocal correlation energy. This was justified because, as a consequence of the rapid electronic motions, the nuclei are immune to the fluctuations of the Coulomb forces, and therefore, the charge distribution must deform to produce the required forces on the nuclei by classical Coulomb interactions alone. These forces can be calculated, as shown by Hellman[73] and Feynman,[74] by exploiting the stationary property of the energy with respect to variations in the wave function. Our calculations with the vdW-DF and revPBE+LDA functional allow for a similar treatment of the molecules considered in this work. In the adsorption of isooctane, the nonlocal interactions cause a nonisotropic accumulation of charge in the region between the molecule and the slab (Figure 4). The electrostatic forces arising from this charge redistribution are responsible for the adsorption of isooctane.

In addition to dispersion forces, a weak electrostatic interaction contributes to the adsorption of ethanol on Fe(100). According to a comparison between the results obtained with revPBE+LDA and vdW-DF, the amount of charge that accumulates as a consequence of nonlocal interactions is larger and extends over a wider area (Figure 5). However, in contrast to the case of isooctane, in this system the charge redistributions predicted by these functionals differ considerably. According to the results obtained with the revPBE+LDA functional, in the absence of nonlocal correlations the Pauli repulsion produces a region of deficit of charge between the surface and the molecule and displaces the electronic density above the molecule. This considerable redistribution of charge (Figure 5a) and the small adsorption energy (Table 1) calculated with this functional suggest that it is overly repulsive in this system. A similar behavior has been observed before for the revPBE functional (exchange and correlation) in molecules at short separations.[75] A calculation with the revPBE functional predicts a larger absorption energy (Figure 8), but it is nevertheless less than half of the one calculated with the PBE functional. Evidently, this system is heavily influenced by the description of the exchange energy, since the only difference between these two functional lies in the exchange term. The charge redistribution calculated with the PBE functional (Figure 6a) resembles the one calculated with vdW functionals (Figures 5b and 6b). Since PBE does not include the effects of nonlocal correlations, this last result indicates that an important contribution to the binding between the ethanol and the Fe(100) surface cannot be attributed to the dispersion forces alone. On the basis of previously estimated adsorption energies, it has been proposed that a weak chemisorption is involved in the binding mechanism.[56] Our calculations, however, do not show a charge transfer large enough to consider the formation of an ionic bond, and the analyses of the electron localization function (ELF)[76,77] do not point to the existence of localized electrons forming a bond between the molecule and the metallic slab (Figure 7). Furthermore, it has been shown that the oxygen states are almost unperturbed,[56] since they do not exhibit a strong hybridization with the d states of the binding Fe atom. The optB86b-vdW functional predicts a lesser amount of charge between the ethanol and the iron slab than the one calculated with PBE, particularly between the hydroxyl group and the closest iron atom to it (Figure 6c). Since this effect is accompanied by an increase in the adsorption energy, a weak electrostatic interaction together with the dispersion forces may constitute the main mechanisms contributing to the adsorption. Nevertheless, this comparison does not allow separation of these contributions, since these functionals differ in the exchange term.

Figure 5

Charge density difference (ρdifff) of ethanol adsorbed on the bcc Fe(100) surface at d = 2.00 Å. The charge density difference is defined and plotted in an analogous way to Figure 4.

Figure 6

Charge density difference (ρdiff) of ethanol adsorbed on the bcc Fe(100) surface at d = 2.00 Å. The charge density difference is defined and plotted in an analogous way to Figure 4.

Figure 7

Isosurface at ELF = 0.50 of ethanol adsorbed on Fe(100) at d = 2.00 Å calculated with the optB86b-vdW functional. Electrons “outside” this isosurface are expected to be delocalized.

The difference in the adsorption process of isooctane and ethanol on a bcc Fe(100) surface can be understood in terms of the polarizability of the molecules and the charge density distribution around the functional group. Clearly, isooctane is expected to be more polarizable than ethanol because of its larger molecular size. Additionally, it is well-known that alkanes are among the most polarizable molecules.[78] Since dispersion forces arise from the formation of instant multipoles, their strength is directly related to the polarizability. Moreover, the absence of a functional group in isooctane leaves dispersion forces as the only possible binding mechanism. On the other hand, ethanol has a permanent dipole, product of the difference in electronegativity between the hydrogen and the oxygen atom in the hydroxyl group, and nonbonding electrons in the oxygen. During the adsorption process, not only the dispersion forces contribute to the binding but there is also a weak electrostatic interaction between the hydroxyl group and the iron slab, where the charge redistribution leads to the formation of multipoles on the top layer (Figures 5 and 6). Consequently, the iron surface binds ethanol more strongly than isooctane even though the contribution of the dispersion forces is expected to be weaker than in isooctane.

Summary and Outlook

The vdW forces are essential for the adsorption of isooctane and ethanol on a bcc Fe(100) surface. As product of these long-range interactions, the nonlocal correlation leads to an increase in the adsorption energies and a reduction of the equilibrium distances. Nevertheless, they do not influence the spatial configuration of the adsorbed molecules. Their consideration in the calculation of the electronic density leads to a more attractive potential in the space between the molecule and the slab. Consequently, less charge is expelled from this region.

Our calculations are a first approach toward a more rigorous treatment of vdW interactions in complex systems, and they will contribute, after experimental validation, to the development and improvement of vdW functionals which are independent of external input parameters. We showed the effects of nonlocal interactions on the electronic density and on the adsorption mechanism. Isooctane binds to the Fe(100) surface via dispersion forces, while in ethanol, in addition to the dispersion forces, a weak electrostatic interaction between the hydroxyl group and the iron surface contributes to the binding.

We anticipate that with the continuous increase in computing power, future calculations that consider many-body effects in combination with the work here presented will clarify the relevance of these interactions and contribute to a more accurate analysis of the vdW forces in systems of industrial interest. In particular, the obtained results will aid the fitting of accurate surface–lubricant interaction potentials required for classical MD simulations including numerous organic molecules. These interface potentials usually constitute a considerable uncertainty in any such MD simulation, as they are rarely assumed more advanced than a Lennard-Jones potential parametrized according to desorption data. A detailed knowledge of the interaction energy between molecules and metal surface can therefore greatly boost the precision of large-scale atomistic studies of the thermal, mechanical, and structural stability of molecular surface films.