Molecule-decorated rutile-type ZnF2(110): A periodic DFT study

Author’s Accepted Manuscript CO-Decorated Rutile-Type ZnF 2(110): A Periodic DFT Study Ali Abbaspour Tamijani, Elham Ebrahimiaqda

www.elsevier.com

PII: DOI: Reference:

S0039-6028(16)30665-3 http://dx.doi.org/10.1016/j.susc.2017.03.002 SUSC20997

To appear in: Surface Science Received date: 16 November 2016 Revised date: 27 January 2017 Accepted date: 4 March 2017 Cite this article as: Ali Abbaspour Tamijani and Elham Ebrahimiaqda, CODecorated Rutile-Type ZnF 2(110): A Periodic DFT Study, Surface Science, http://dx.doi.org/10.1016/j.susc.2017.03.002 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

CO-Decorated Rutile-Type ZnF2(110): A Periodic DFT Study Ali Abbaspour Tamijania*, Elham Ebrahimiaqdab a

Department of Chemistry, Wake Forest University, Winston-Salem, NC 27109 USA Department of Chemical Engineering, University of Arizona, Tucson, AZ 85721 USA

b

*

Corresponding author: [email protected]

Abstract Weak binding of small molecules onto surfaces is a powerful tool whereby interfacial phenomena can be studied in atomistic level. It is well-recognized that physisorption is a precursor to chemisorption and the subsequent heterogeneous catalysis. Although at least in part overlooked, vdW-driven sorption of particles on rutile-like substrates is of potential value due to the wide variety of the applications it serves. Probing the acidity of rutile-structured adsorbents by means of carbon monoxide adsorption is the quintessential case of long-rangedominated adsorption on such materials. Monomer, half-layer and monolayer physisorption of gaseous CO at (110) facet of rutile-like ZnF2 was investigated through dispersion free and dispersion-corrected DFT. Commensurate with the pertinent literature, the upright C-ligated monomer when bound to unsaturated Zn, was found to be energetically favored over the Oligated conformer. The PBE and optB88-vdW calculated stretching frequencies for freestanding and adsorbed CO were in excellent agreement with their experimentally observed counterparts. Traditional vdW-DF and its successor, vdW-DF2, predict surface energies that match well with B3LYP result. As well, outstanding consistency was found between the vdWDF and vdW-DF2 computed binding energies and their highly accurate LMP2 equivalent. Whilst being computationally far more efficient, vdW-DF and vdW-DF2 attained interaction energies were closer to the LMP2 benchmark data than the previously reported B3LYP adsorption energy. Keywords: Atomic displacements, adsorption energy, dispersion-corrected DFT, CO, rutile, ZnF2.

1

I. Introduction -Background: Substances interact with their surroundings through surfaces. To gain a deeper insight into materials and their functionalities, a thorough investigation of the microscopic structure of surfaces is required. For a long time, the scientific community’s knowledge of the structural details of the rutile-like binary compounds was confined to oxides. Aside from rutile-like oxides, their fluoride analogues are worthy of attention due to the variety of their remarkable applications. For instance, metal fluorides provide interesting properties such as significant acidity [1] and wide band gap [2]. As a transition metal fluoride, ZnF2 is regarded as a well-known optical compound. It has driven much interest as a potent electroluminescent material [3]. It is also well-documented that zinc fluoride acts as an attractive subject of study for vacuum ultraviolet photo-electron spectroscopy due to its wide band gap [2]. This material has recently drawn considerable attention due to the discovery of new sol-gel synthetic strategies that can beget nanostructured metal fluorides with large surface areas [1]. Metal fluorides are frequently utilized in anti-reflection coatings and in ceramics industry [3]. Additionally, zinc fluoride’s diverse dielectric properties upon doping of Li has made it the target of primary significance in optical device manufacturing [3]. Although overshadowed by the well-deserved fame of the iconic rutile TiO2 and similar transition metal oxides, quasi-rutile zinc fluoride can also function as heterogeneous catalysts [3]. Given the recent developments in the synthetic strategy of highly catalytic magnesium and zinc fluorides [1], the surfaces of these materials have given rise to a pronounced interest amongst tribology scientists [4-6]. The synthetic rutile MgF2 has exhibited both outstanding catalytic activity and selectivity. Based on the similarities that coexist between the electronic structures and crystalline behavior of MgF2 and ZnF2, the latter has also become the subject of ongoing research [1]. Given the aforementioned motivations, the electronic structure description of ZnF2 surface appears to suffer from a gap due to insufficient literature resources. Thus far, there has been only one study that has revealed the surface energetics and structural details of several lowindex facets of rutile zinc fluoride [3]. The abovementioned study has exploited B3LYP along with extended Gaussian-type basis sets implemented in the CRYSTAL09 program. As with any other method, a hybrid scheme has both advantages and disadvantages of its own [7]. When treating systems with predominantly non-covalent interactions, it is a hurdle to unequivocally assess the quality of adsorption energies computed via B3LYP. In some instances, B3LYP computed binding energies have been improved with the aid of common dispersion correction schemes such as D3. Examples are adsorption energies calculated for H2-bound graphene, COcovered hydroxyapatite nanoparticle surface and alkylamines at H-mordenite [8-10]. On the other hand, non-dispersion-corrected B3LYP has shown promise for certain adsorption cases. For instance, the adsorption of CO atop the surface of Co3O4 was simulated using this approach alongside effective core potential (ECP) basis sets and the results were justifiable [11]. The same level of theory and basis functions were employed to model the sorption of water atop low-index Pt clusters and chemisorption energy predictions were satisfactory [12]. Hence, we were incentivized to perform periodic relaxation calculations on ZnF2 using dispersion-free and dispersion-corrected approaches other than B3LYP. This will help understand the effect of van der Waals forces on surface structure and energy. Given the better computational efficiency of a non-hybrid exchange-correlation functional, this will also allow for the evaluation of performance of such approaches compared with a hybrid one. The most stable crystalline form of ZnF2 transpires in the form of rutile [3]. As expected, this polymorph exists under ambient conditions. Amongst low-index planes of rutile-like ZnF2, (110) is found to be the most energetically preferred [3]. The Kaawar and Paulus study 2

has reported a surface energy of 0.40 Jm-2 for the (110) plane of this material [3]. This value is closer to the lower limit of the surface energies reported for other (110) facets of rutile-type minerals [13,14]. This is an indication that the energy penalty to create this surface is relatively low. This is attributed to the weak nature of the Zn-F bond. This effect itself is by virtue of the electronic structure of ZnF2 being dominantly of covalent rather than ionic nature [15]. Another reason for the comparatively low surface energy of zinc fluoride is the flexibility of the axial Zn-F bond. The literature is contradictory on whether the apical Zn-F bond favors elongation or compression. Nevertheless, the small energy barrier for such a structural change can facilitate the bond-clipping process, resulting in a more stable surface [16]. Decorating the surfaces of crystalline solids with atoms and small molecules is identified as a useful means whereby properties such as substrate acidity can be explored [1,17]. Among the molecules used to probe the surfaces, CO, N2 and NH3 are of paramount value due to their increasing popularity in such techniques [6]. In addition, it is well-recognized that adsorption is commonly a precursor to subsequent phenomena such as heterogeneous catalysis. Amongst the forces that drive adsorption phenomena, dispersion effects are particularly noteworthy. These weak effects occur in nature frequently. In spite of their prevalence, these forces are difficult to treat. Physisorption is mainly governed by classical electrostatic and/or dispersion interactions [5]. This type of adsorption has found many applications in tribology such as monitoring the growth of adlayers onto the surfaces of solid substrates [18]. Examining the literature, countless prototypical examples of physisorption such as weak binding of small organic molecules to insulating ionic crystals have been theoretically studied at length in the past [19-21]. The investigation of the adsorption of CO on the surface of β-AlF3 was carried out by Krahl and co-workers [22]. The work has put forth metal fluorides as a new class of heterogeneous catalytic materials. In order to study the Lewis acidity strength of aluminum fluoride, its CO-physisorbed surfaces were spectroscopically scrutinized by IR instrumentation. The resulting spectra were examined and the adsorbate-substrate vibrational modes were identified. Another article, which also has inspired the current work, contains the results of a set of experiments conducted by Guo et al. [1]. According to Guo and co-workers, small organic molecules such as pyrrole, chloroform and carbon monoxide are to be used as probe molecules to characterize ZnF2. First, they synthesized the zinc fluoride by a fluorolytic sol-gel technique. Then, they conducted highly sensitive adsorption experiments. Afterwards, Guo and colleagues used in situ FTIR to explore the surface of the yielding material for its acidic and basic properties. They have assigned the peaks of the resulting spectrum to vibrational motions of the adsorbate-adsorbent couple. The abovementioned experimental results have been recently confirmed through theoretical investigations done by Kaawar and co-workers [17]. They modelled the adsorption of CO on top of rutile ZnF2 surfaces with several lowindex structures. B3LYP results were found to be acceptable when compared with benchmarking LMP2 binding energy values. Buchholz et al. has outlined that infrared reflection absorption spectroscopy can characterize ad-CO probe particles on single-crystal surfaces of ZnO applying p-polarized and s-polarized radiations. To this end, close to upright geometry was found for adsorbed CO atop the Zn2+ of the zinc oxide substrate. Coverage-dependent vibrational frequency change was discovered for ad-CO stretching mode. It was concluded that upon increasing the surface coverage, the vibrational motion frequencies become red-shifted [23]. -Methodology Density Functional Theory (DFT) has put forward a viable path to successfully describe the electronic structure of systems of many electrons [24]. In recent years, this method has become 3

the workhorse for modelling condensed matter systems [25]. Despite its sophistication, traditional DFT does not include van der Waals effects. This can be problematic when treating systems of a dispersion-dominated nature. To circumvent this predicament, several approaches have been introduced. a. DFT-D3: According to this method, the dispersion term is additive and the total energy can be expressed as [26]: EDFT-D3=EKS-DFT+Edisp

(1)

In equation (1), EDFT-D3 refers to the total dispersion-corrected energy, EKS-DFT is the KohnSham orbitals energy value obtained by DFT, and Edisp denotes the dispersion energy. Edisp is assumed to be a sum over two- and three-body terms: Edisp=E(2)+E(3). The inclusion of a damping function is of absolute necessity for this scheme. The damping function can be of the form proposed by Becke and Johnson (BJ) [27]. It can also be the zero-damping function. Several DFT (e.g. PBE [24]) and hybrid (e.g. PBE0 [28]) methods are corrected for weak interactions using this approach. It has been shown that the effect of three-body interactions can be up to 25% of the total energy for some systems [29]. These effects are therefore of particular value for condensed materials and sorption cases. b. vdW-DF family: Originally proposed by Dion et al. in 2004 [30], it treats the exchangecorrelation contribution as:

E xc  EGGA  EcLDA  Ecnl  x

(2)

In equation (2), denotes the exchange-correlation part, while is the exchange-only term computed via the generalized gradient approximation. refers to correlation term, which is calculated through the local density approximation. is the non-local correlation effect, given by:

Ecnl 

1 3 3 d rd rn  r  Φ0  r, r  n  r  , 2

(3)

where n  r  and n  r  denote the electron densities at points r and r , respectively. The kernel

Φ0  r, r  is a universal function of the n  r  ’s and their gradients. The local correlation term is chosen as LDA. The exchange term in the non-local method is in the form of GGA. Several variants of vdW-DF methods exist and are actively in use for different systems. For example, vdW-DF which is the initially proposed non-local approximation, uses revised PBE (revPBE) to account for the exchange energies [31]. This exchange functional differs from its predecessor (i.e. PBE) in only one variable of its exchange portion, viz., κ, which has been altered from 0.804 in PBE to 1.245 in revPBE [32]. With modernized implementation of vdWDF and its family, the level of computational effort involved in these calculations has become more practicable. vdW-DF2 uses an improved exchange functional (rPW86) compared to vdW-DF [33]. The vdW-DF2 has fared well in several instances. As an example, the adsorption of n-butane on Cu(100), Cu(111) and Pt(111) can be mentioned for which adsorption energies predicted by vdW-DF2 are in excellent agreement with experiments [34]. Ma et al. [35] has tested the performance of vdW-DF2 for CO-adsorbed transition metal surfaces for twenty-five systems comprised seven transition metal elements. It was deduced that in its present implemented algorithmic form, vdW-DF2 makes satisfactory predictions of adsorption energetics and vibrational frequencies for the systems of interest. Due to its notably good results for CO 4

adsorption, the vdW-DF2 method has also been used in the current work. Other implemented versions are optPBE-vdW, optB88-vdW and optB86b-vdW, with the second and the third performing similar to each other for tested systems [36 and references therein]. Each non-local variant is particularized to subdue a specific deficiency. This approach exploits an optimized version of the well-known B88 exchange functional. Its adsorption energies have shown consistency with temperature-programmed desorption (TPD) experimental data for benzene adsorbed onto (111) surfaces of transition metals [37]. Aside from all these promising results, this method, has been previously applied to rutile TiO2(110) and has performed reasonably. In the study done by Tillotson et al. [38], optB88-vdW surface energy and surface-induced atomic displacements of r-TiO2(110) have been calculated and validated when compared to experiments. As well, adsorption energies for methane, formic acid and glycine (all being small organic molecules as is CO) have been calculated at this interface. These encouraging results led us to believe that optB88-vdW is also a good choice for examining the (110) face of r-ZnF2, a crystalline surface of geometrical parameters comparable to those of r-TiO2(110) one. The present paper focuses on the investigation of the contribution of long-range interactions to the adsorption energies of CO molecules on the surface of rutile-type zinc fluoride. We have examined the effect of lateral interactions on the overall adsorption energetics. The vibrational frequencies of ad-CO have been calculated and compared to the experiments. The framework of the current article consists of the following. Firstly, the structural details of the relaxed (110) surface of ZnF2 is reported for GGA-PBE and optB88vdW methods. Secondly, the adsorption geometries and energies for low-coverage COadsorbed surface will be discussed. Finally, the effect of lateral interactions for half- and mono-layer adsorptions are studied as well as the impact of adsorbate-induced surface relaxations.

II. Computational Details The parallelized version of the commercial Vienna Ab initio Simulation Package (VASP) program was used [39,40]. A Monkhorst-Pack [41] 2×2×1 grid of k-points was chosen to sample the BZ. Kohn-Sham basis sets were expanded in terms of plane-waves of up to 700 eV of kinetic energy cutoff. A 4×2×4 supercell consisting of 192 atoms was chosen. The size of this semi square-shaped surface ensures the convergence of the surface energies obtained for the material of interest. The experimental lattice constants listed in Table 1 were used to create the supercell [42]. A vacuum gap of at least 17 Å of thickness was inserted above the slabs to avoid the surface interaction with its periodic image along the c crystallographic axis. Forces were corrected by adding the ADDGRID=.TRUE. command line. Pseudopotentials used for Zn, F, C and O comprised 12, 7, 4 and 6 valence electrons, respectively. A 0.05 eV wide Gaussian smearing was introduced to account for partial occupancies. All of the surface atoms were allowed to move during the relaxation calculations. However, the slab shape and volume were preserved. ZnF2 surface contains cations and anions that have full valence orbitals (either p or d). Zn possesses a +2 oxidation state in ZnF2 and its fully occupied d orbitals show diamagnetic properties. This has been corroborated by 0.00 magnetic moments calculated for ZnF2 using both LDA and PBE methods [43]. Thus, no spin polarization effect was considered in the calculations. A choice which boosts computational tractability. GGA method proposed by Perdew-Burke-Ernzerhof (PBE) was also exploited [44].

5

Table 1. Experimental values of structural parameters of rutile-type ZnF2 [42].

a[Å]

c[Å]

u

4.703

3.133

0.305

To include dispersion effects in the surface structure, optB88-vdW was employed to relax the slab. This method has been previously applied to rutile TiO2(110) and has shown promise in prediction of both surface energy and atomic displacements [38]. Zinc fluoride coordinates were frozen at the optB88-vdW level for all dispersion-corrected methods. For the sake of consistency, the geometry of the CO molecule was also kept at the same level throughout the body of this work. Minimizing the slab and adsorbate structure with one method and performing the binding energy calculations with another, is not an unprecedented approach. As two examples, Zhao et al. [45] and Bajdich and co-workers’ [46] studies can be mentioned. The former group has investigated the co-adsorption of CO and H2 on Ru(0001) surface and the synergic effects between the two adsorbed species. RPBE calculations were done on a PBE-relaxed slab. The latter work however, is a systematic study on the binding of small particles, namely, CO, NO, CH4 and H2O onto MgO, CaO, SrO and BaO using RPA at the PBE-optimized geometry. In the absence of experimentally measured atomic displacements of ZnF2(110), no choice was left but to utilize a method (optB88-vdW) which has been experimentally validated for a similar surface structure (r-TiO2(110)). For PBE calculations, the coordinates were frozen at the PBE level of theory. PBE and dispersion-corrected results were also used to compute surface energies. To do so, the following formula was applied: (4) In equation (4), , , and are surface energy, slab energy, conventional unit cell energy and surface area of the chosen slab, respectively. is the number of unit cells in the slab. CO molecule(s) was (were) added to one side of the slab. Although the defect-free (110) surface is non-polar, placing the ad-molecule on one face of the slab can result in artificial dipole moments. To eliminate the errors caused by mono-, di-, and quardru-pole moment formation, the Makov and Payne correction was applied [47]. IR frequencies were computed through a harmonic approach. For the adsorbed CO, only vibrations in the direction normal to the plane were considered. One way to quantify the likelihood of the adsorption process to take place is the adsorption energy. All the adsorption energies were computed via:

Eads 

Eads nCO   Eslab  nECO   n

(5)

In equation (5), E ads is the adsorption energy. n is the number of CO’s placed on the slab. Eslab and E CO represent bare slab and isolated CO energies, respectively. The interaction potentials were obtained by computing single-point energies within 2.00-10.00 Å. The potential energy curves were smoothened with a typical physical adsorption potential energy function. Proposed and tested for the sorption of noble gas atoms onto solids [36,48], the function includes a short-range Pauli repulsion term and a long-range asymptotical contribution. The overall form of the function is: ( )

( )

( ),

(6)

where ( ) is the potential energy as a function of the adsorbate-adsorbent distance, z. The term ( ) is the short-range repulsive part. ( ) on the other hand, denotes the long6

range zone of the potential energy curve. The short-range portion of the interaction potential originates from the overlap between the sorbent conduction band and the valence electrons of the adsorbed species [49]. In many cases, the application of an exponential expression of the ( ), with form ( ) and being adjustable parameters, has led to favorable ( ) signifies the long-range interaction end-results in describing the repulsive potential. terms(s). It has been demonstrated that the weak attractive zone of the gas adsorption energy curve can be modelled by and coefficients. The former is shown to depend upon the adsorbate polarizability and the surface work function. The functional behavior of the latter however, is unknown and this term is parametrized using the ab initio potentials. The longrange attractive function was used alongside a damping factor that accounts for the effect of the long-range decay of the potential energy. The following section contains a summary of the gas phase calculations for isolated CO as well as surface energies and atomic displacement results for relaxation calculations using PBE and optB88-vdW. It also presents computed binding energies and adsorption heights for lowto mono-layer coverages using a multitude of dispersion-corrected DFT methods. Together with a discussion and results validation, the predicted values are evaluated by making comparisons to literature data for similar systems.

III. Results and discussion -Freestanding CO Polar heteronuclear diatomic CO has an experimental bond length of 1.1283Å [50] and a reported stretching frequency of 2100cm-1 [51]. The geometry of gaseous CO molecule was optimized in a box using both dispersion-free (PBE) and dispersion-corrected (optB88-vdW) exchange-correlation functionals. The results are presented in Table 2. Table 2. Bond length (Å) and zero-point vibrational frequency (cm-1) of gas-phase CO: theoretical vs. experimental. PBE

optB88-vdW

Exp.

C≡O

1.1355

1.1325

1.1283 [50]

υCO (cm-1)

2125

2115

2100 [51]

While all of the DFT-based values are slightly overestimated, acceptable agreement exists between the experimental and theoretical results. Including the effect of weak interactions in optB88-vdW betters the agreement between the experimentally observed and theoretically predicted quantities by a minute amount. It is expected that GGA-PBE method should underestimate the bond length by 2-3% as this has been observed for CO and other diatomic molecules in the past [52]. Nevertheless, the dispersion-less PBE computed equilibrium separation of CO is in excellent agreement with the experimental value. Calculated vibrational frequencies are all consistent with the experimental data. Other works have reported a vibrational frequency of 2146 cm-1 using a GGA-PBE approach [53] and a value of 2191 cm-1 using the GGA+U functional [54]. As such, the 1-3% difference between our PBE computed value and the literature PBE-based data is certainly justifiable. Compared to the experimental value for gas-phase CO, the optB88-vdW method makes the least erroneous prediction of all (with an absolute relative error of ~0.7%). Our optB88-vdW bond length value of 1.1325Å is only ~0.41% larger the almost error-free value of 1.1279 Å, obtained by Dawes and coworkers using highly precise CCSD(T)-F12b(AE)/CVQZ-F12 method [55]. Kaawar et al. has reported a vibrational frequency of 2206 cm-1 using B3LYP [17]. This value resides within ~3% of the experimental data referenced therein (i.e. 2141 cm-1). Our PBE and optB88-vdW 7

theoretical frequencies deviate from their benchmark value by ~0.75% and ~1.21%, respectively. -Bare ZnF2 Out of all the low-index rutile surfaces, the (110) slab is shown to be the most stable cut. This has been confirmed repeatedly through calculations and experiments for numerous materials exhibiting rutile crystallinity such as TiO2 [56], SnO2 [14], MnO2 [57], MgF2 [58] and ZnF2 [3]. The (110) surface contains two types of cations (i.e. Zn2+), namely, five- and six-fold. It also has two types of anions (i.e. F-), namely, two- and three-fold. The three-fold anion is also referred to as in-plane. Also, the two-fold anion is sometimes labeled as bridging or apex. The two-fold atom is protruding from the surface. Figure 1 is a typical (110) plane of rutile ZnF2. No low-energy electron diffraction (LEED), surface X-ray diffraction (SXRD) information or computational benchmark value is available on the atomic coordinates of the (110) surface of rutile ZnF2. Nonetheless, this data is available for similar systems such as TiO2 [38]. Closely examining these experiments conducted for similar surfaces, it can be inferred that the pattern with which the surface atomic shifts occur are the same for all (110) rutile-type planes [38,56,59]. Compared to the bulk-terminated structure, the under-coordinated cation moves farther into the surface. This downwards motion creates a cave-shaped hollow spot that functions as an acidic site. The four-fold cation moves outwards to adopt a distorted tetrahedral geometry. The remainder of the surface atoms also undergo an upwards motion upon relaxation. As the atomic displacements mainly happen perpendicular to the surface, lateral relaxations are negligible. Accordingly, this work is only concerned with atomic shifts along the c crystallographic axis. Second and third layers of rutile-like surfaces experience displacements too. Although, due to the presence of missing bonds in the experimental measurements, the amount with which these changes take place are different than theoretical values as the latter is mainly assuming perfect crystal structures, an approximation that can be misleading in comparison with the reality. However, both theory and experiments predict that the amounts of relaxations erode with depth. In order to assess the role of k-point grid density in the relaxation calculations, both a Γcentered single k-point and a 2×2×1 Monkhorst-Pack mesh were employed to relax the slab configuration. We concluded that the density of k-points does not have a pronounced effect on the displacements. Table 3 summarizes the obtained values for a 2×2×1 grid. In this Table, a negative displacement denotes an inwards movement whereas a positive one indicates an outwards shift (towards the direction of the vacuum gap). Zn5f, Zn6f, F2f and F3f are coordinatively unsaturated Zn, second layer six-coordinated Zn, bridging and in-plane F atoms, respectively. Table 3. Atomic displacements of the two topmost layers of the ZnF2(110) along the c direction [Å]. Zn5f

Zn6f

F2f

F3f

PBE

-0.03

0.40

0.27

0.27

optB88-vdW

-0.08

0.33

0.19

0.18

As far as the geometrical details go, the surface of zinc fluoride experiences the same qualitative level of change as the rest of the rutile-like (110) planes. As expected, coordinatively unsaturated Zn (Zn5f) moves further into the bulk and forms a nested area that 8

resembles a square pyramidal structure with its basal plane facing upwards (towards the vacuum). Irrespective of the choice of the exchange-correlation functional form, the general trends of the atomic displacements follow the same pattern as above. The amount of displacement however, is sensitive to the type of the functional employed. Our results indicate that on accounting for the effect of weak interactions, under-coordinated zinc undergoes a greater displacement. The amount of displacement for five-fold zinc when computed through optB88-vdW is almost tripled compared to the PBE prediction. Regardless of the method, twoand three-fold fluorides experience very similar shifts along the direction normal to the (110) plane. Overall, the values of displacements agree well with those observed for other rutile-like (110) facets. Zn-F bond distance is close to the sum of cation and anion radii (i.e. 0.74Å and 1.33Å for Zn and F-, respectively [60,61]). The axial Zn5f-F bond for the bulk-terminated version of the slab is 2.03 Å, which decreased to 2.00 Å for the PBE-computed relaxed version. As a percentage, the bond length reduction (~1.5%) is much less than in oxides such as SiO2 (6%) [62]. As such, the change in displacement is much smaller than that observed for rutile oxides. The axial distance reduces to 1.98 Å upon relaxing with the optB88-vdW functional. The equatorial Zn5f-F3f bond for the bulk is 2.03 Å, whereas for the PBE-calculated relaxed slab, it is 2.02 Å, as is also the case when optB88-vdW is employed. This is consistent with previous observations for ZnF2 in its bulk form [43]. As opposed to oxides [62], the axial Zn-F bond remains longer than the equatorial Zn-F bond. Zn6f-F2f distance is 2.03 Å before and after relaxation, respectively when computed via PBE. The relaxed distance calculated using optB88-vdW is 1.95 Å. The stiffening of this bond is a consequence of non-equidistant motions of Zn6f and F2f upon relaxation. As expected, the distance between two neighboring Zn5f atoms (ions) is 3.13 Å for both F-terminated bulk and PBE-optimized slab as well as its optB88-vdW-relaxed counterpart. This distance is slightly longer than that of TiO2 (2.96 Å) [63]. Since this distance is known to be a determining variable in some of the TiO2’s catalytic activities (such as N2O2 dissociation into 2NO’s [64]), ZnF2 can also be tested for such functionalities as its structure shows similarity to the titanium dioxide surface. The bond angle Zn6f-F2f-Zn6f is 100.7º and 106.7º before and after geometry optimization with PBE, respectively. The optB88-vdW functional predicts a very slight enlargement of this angle to 106.9º. Several literature values were found for F3f-Zn5f-F3f including, but not limited to, 79º [16] and 80.51º [43]. Both of these values confirm that the bulk structure is a distorted rutile. Dispersion-less displacement of five-fold Zn is underestimated as for similar systems, this value is ~0.10-0.20 Å (as is for TiO2 [38], SnO2 [59] and FeF2 [65]). This distorted octahedral environment concomitantly facilitates the (110) surface formation. As mentioned before, the F-Zn-F angle between the axial and the equatorial in the bulk is 80.51º [43]. Our PBE-relaxed counterpart of this angle is 80.39º, in satisfying agreement with the experiments. 2+

In the absence of a consensus about the exact amounts of atomic displacements, no final decision can be made as to whether or not any of the selected methods is the ‘best’ option for these systems. The scarcity of literature data also leads to incapability of assessing the impact of dispersion interactions on the structure as a whole. It is essential to keep in mind that even with experimental data available, the difference between the real conditions and approximations made to perform simulations will remain an open issue. Our PBE-estimated surface energy is 0.30 Jm-2 which is distinct from its B3LYP computed counterpart. The disparity is clearly a consequence of the differences in treatment of the exchange-correlation density functional with the two methods. It is recognized that surface energies in a periodic approach are reliant on the type of the exchange-correlation functional as well as the slab thickness [13,66]. The low PBE surface energy value of 0.31 Jm-2 for TiO2(110) supports the notion that this method predicts lower surface energies for similar crystals [13]. Our optB88-vdW computed value is 0.62 Jm-2. The optB88-vdW calculated 9

surface energy differs from that using B3LYP in the virtue of the same two reasons described above. It is essential to emphasize that the addition of the dispersion effects (regardless of the dispersion inclusion approach) almost always increases the surface and binding energies compared to dispersion-less PBE method. The optB88-vdW computed surface energy is larger than that of the dispersion-free hybrid method. The difference between PBE and optB88-vdW is profound. Table 4 itemizes all of the computed surface energies. Table 4. ZnF2 surface energies (Jm-2) obtained via different methods.

PBE

vdW-DF

vdW-DF2

optB88-vdW

B3LYP [3]

0.30

0.41

0.42

0.61

0.41

Both vdW-DF and vdW-DF2 predict surface energies in excellent agreement with the B3LYP-computed value, reported by Kaawar et al. [3]. The optB88-vdW overestimates the energy by ~49% whereas PBE underestimates it by ~37%. Since the only difference between variants of nonlocal approach is the way the exchange terms are treated, it is concluded that the separation between the optB88-vdW surface energy result and the vdW-DF and vdW-DF2 computed values is given rise to by inherent differences in their exchange energy estimations. Due to the lack of experimentally measured surface energy values at low temperatures, neither our results nor the previously reported energy value can be directly examined and compared for their veracity. To accomplish this task, further investigations are required in form of either experiments or high-level ab initio benchmarking calculations (such as LMP2). Until then, the only means whereby we can ratify our results is through literature values for similar systems. However, making comparison to similar systems has its own deficiencies such as not accounting for the effect of possible anomalies that can occur among these systems and uncertainties imposed by other tangential factors. -CO adsorption Physisorption is the result of a delicate balance between repulsive and attractive forces. CO adsorption on crystal facets has long been the subject of numerous experimental and theoretical studies [67-70]. Besides TiO2 [71], this small molecule is proven to be an aid in characterization of several other rutile-structured surfaces [4-6]. The main goal of this segment of the work is to model the adsorption of CO onto rutile-like ZnF2(110) surface. To serve this purpose, computation of the binding energy of CO is imperative. The binding energy indicates the strength of the bond formed between the adsorbate and the surface. It also helps find the most energetically preferred configuration of the ad-species. To avoid the unnecessary computational cost, only the upright geometry of adsorbed CO atop the under-coordinated Zn was considered. It has been well-documented that this conformation is the most stable of all possibilities for CO at TiO2(110) [67], SnO2(110) [72], RuO2(110) [73], MgF2(110) [5,6] and more recently, ZnF2(110) [17]. Accordingly, only two potential isomers are examined hereafter: C- and O-bound. To the best of the authors’ knowledge, there is only one set of robust quantum chemistry-computed adsorption energy results available for the CO-ZnF2(110) binary system. The research conducted by Kaawar and co-workers [17] has outlined the adsorption energy values for 0.5ML and 1.0ML of coverage using both B3LYP and LMP2.

High-level incremental MP2 and CCSD(T) calculations have been performed for COMgF2(110) system and have predicted the adsorption energy to be 373-385meV, depending on the method exploited [5]. MgF2(110) and ZnF2(110) are similar in many aspects such as their electronic structure and crystalline arrangements as well as their ionic radii (0.72Å vs. 0.74Å 10

for six-coordinated Mg2+ vs. six-fold Zn2+, respectively [60,61]). However, the surface of ZnF2 is less acidic than MgF2. This originates from the fact that the electronegativity of Zn (χ=1.65) is higher than that of Mg (χ=1.31). As such, the Zn-F bond is less polar than the Mg-F one. This effect leads to a denser charge distribution around Zn compared to Mg, which itself will make Zn less affinitive towards ad-species. The lesser acidic strength of Zn means a smaller binding energy [1]. This effect has been interpreted as the result of partial covalent nature of the Zn-F interaction [15]. Setting the CO-MgF2(110) adsorption energy value as the threshold, the CO-ZnF2(110) counterpart ought to be less than the energy range mentioned above. The differences must be small however, as the acidities are not expected to be dramatically different. Notwithstanding that there is no experimental evidence that leads to a binding energy, spectroscopic data is in fact reported for the vibrational frequency of CO physisorbed on zinc fluoride surface. Guo and co-workers have outlined a range of 2158-2153 cm-1 for the vibrational motion of CO atop five-fold Zn [1]. This quantity is extracted from the differential IR spectrum presented therein. All the frequencies are greater than the free monomeric CO frequency of 2100 cm-1. The frequency shift is explained by Angel and Schaffer’s work on CO adsorption onto zeolites [74]. The work has put forth the undoubtful existence of a correlation between the frequency of adsorbed CO and the electric field in the vicinity of the surface cationic site. The influence of this electric field on the carbon nucleus of C-ligated CO is proposed to result in the partial removal of the electron cloud from 5σ molecular orbital of CO. Ab initio calculations have demonstrated that the 5σ orbital is of antibonding nature and is primarily situated at the carbon end of the adsorbate [75]. The dislocation of the electron density of the 5σ from its original position towards the substrate metallic ‘hollow’ site, leads to strengthening of the C≡O bond and the subsequent increase of its designated stretching frequency. Consequently, it is plausible to expect the bond length to shorten as CO is adsorbed onto the surface. The increase of the stretching mode frequency of CO is in agreement with the trends observed for CO adsorption on CeO2(111) and anatase TiO2(101). Done by Reimers and coworkers [54] using LDA, the former has reported the vibrational frequency value of adsorbed CO as 2155 cm-1. The latter study on the other hand, has utilized an experimental technique called infrared reflection absorption spectroscopy (IRRAS). This last study revealed a stretching frequency of 2181 cm-1 for CO when adsorbed atop the (101) plane of anatase TiO2 [76]. Conducted by Setvin and colleagues, the same work has highlighted an experimental adsorption energy value of 0.37±0.03 eV for CO-TiO2(101) [76]. Due to the inherent resemblances that exist between transition metal oxides and fluorides, both the blue-shift of the CO vibrational frequency and the amount of the adsorption energy are similar to those observed for CO on top of ZnF2(110). Putting stronger emphasis on the commonality of this behavior, the blue-shifted frequency for the adsorbed CO is also observed for CO-Fe3O4(111) and ZnO. As Lemire et al. [77] has suggested, the adsorbate’s vibrational motion for COFe3O4(111) is accompanied by a frequency of 2115-2140 cm-1, 2080 cm-1 and 2207 cm-1 at temperatures 110K, 180K and 230K, respectively. Additionally, de Lima et al. GGA-PW91 calculations have conveyed that CO vibrational frequency changes to 2183 cm-1 upon adsorption at ZnO(10 ̅ 0) surface [78]. Clearly, all these frequencies are above the 2100 cm-1 value, measured for free CO molecule. The CO vibrational motion on MgF2(110) corresponds to a frequency blue-shift of 57 cm-1 [6]. The CO adsorption onto ZnF2 surface is accompanied by a blue-shift of the frequency within 53-58 cm-1, a slightly smaller shift on average compared to that of MgF2 one [17]. At 0.1ML coverage limit, the stretching frequency of the physically adsorbed CO was predicted to be 2132 cm-1and 2135 cm-1 using PBE and optB88-vdW functionals, respectively. Although both values match the experimental frequency range of 2158-2153 cm-1, the latter affords a 11

better agreement with the information given by in situ IR spectrum. Although inconsiderable, incorporating the weak effects modifies the agreement with the experiments. The CO experimental vibrational frequency shows ~2.52-2.76% increase upon adsorption [1]. Our PBE and optB88-vdW frequencies exhibit a 0.33% and a 0.95% of increase, respectively. The current work has employed several DFT-based approaches to predict both adsorption energy and equilibrium distance between the CO and surface. Table 5 lists all the obtained binding energies for 0.1 ML end-on and head-on isomers. The optimal CO-surface interatomic separation is defined as the distance between the C and Zn5f. In this Table, BJ stands for a Becke-Johnson damped calculation while no mention of the damping function type means a zero-damped calculation. Table 5. Adsorption energies (meV) and, in brackets, equilibrium distances (Å) at 0.1ML for C- and O-ended isomers.

C-bound O-bound

PBE 247.2[2.52] 71.65[2.86]

optB88-vdW 447.0[2.37] 276.0[2.57]

vdW-DF2 364.8[2.63] 233.4[2.68]

vdW-DF 363.0[2.69] 256.4[2.75]

RPBE-D3 377.6[2.65] 225.3[2.77]

RPBE-D3(BJ) 422.5[2.62] 245.5[2.76]

It is demonstrated that regardless of the picked method, the C-terminated conformer is always more stable than the O-terminated. It is inferred that PBE and optB88-vdW approaches underand overestimate the binding energy, respectively. Irrespective of the chosen dispersion correction scheme, inclusion of the weak effects significantly enhances the binding strength. The largest difference between the C- and O-terminated binding energies was observed for the dispersion-free PBE functional. The PBE evaluated adsorption energy is more than tripled when the OC is flipped to CO. The BJ-damped RPBE-D3 method overestimates the binding strength similar to optB88-vdW. However, the zero-damped version is more realistic and close to vdW-DF2 and its predecessor, vdW-DF. The exchange functionals employed in vdW-DF and RPBE-D3 implementations behave close to each other. The resemblance between the binding energies predicted using these two methods dictates that their correlation terms should, for this specific case, also behave similarly. This is particularly interesting since two fundamentally distinct approaches perform at the same level of accuracy regardless of the way the correlation energy term is treated. The energy separation between the C- and the O-bound conformers is observed to be very close for the vdW-DF, vdW-DF2 and RPBE-D3 approaches, where the adsorption of the former is only ~30% stronger than the latter. Nonetheless, all of the binding energy differences are significant for C- and O-terminated adsorbate and the C-bound isomer is always energetically preferred. The RPBE-D3(BJ) adsorption energy shows the biggest change when altering the C-bonded to the O-bonded isomer (177 meV). Upon changing the ligation end of the CO, the binding energy experiences an average reduction of ~152meV. Addition of dispersion effects increases the binding energy for O-terminated isomer as well. Going from PBE to vdW-DF2, the O-appended adsorption energy is more than tripled. The change from vdW-DF2 to optB88-vdW is much smaller. Nonetheless, the binding energy increases by ~18%. The increase is more substantial from PBE to optB88-vdW (~286%). Aside from the PBE calculated binding energy, all other binding energies are rather close for the O-ligated conformer. At low coverage limit, the Cbonded adsorption energies range from 247.2 to 447.0meV. These values are within typical physisorption energy range. An increase is observed for O-ligated adsorption energy predicted going from PBE to optB88-vdW. This variation is a result of the exchange correlation functional choice. Compared to the RPBE-D3, the PBE underestimates the adsorption energy by ~53%. For C-bound isomer, traditional vdW-DF and its amended version, vdW-DF2, generate adsorption energies that are very similar to one another (within ~0.5% of each other). The similarity between vdW-DF and vdW-DF2 results implies that revPBE, as the vdW-DF exchange term, and rPW86, as the vdW-DF2 exchange energy term, perform almost identically for CO-ZnF2(110) system at low surface coverage limit of 0.1ML. As it was 12

discussed earlier, the CO-ZnF2 couple binding energy should be slightly smaller than 373385meV. As a consequence, both of the abovementioned approaches give rise to trustworthy results for the system of inquiry. The other hand, the semi-local RPBE, when applied alongside D3 dispersion correction scheme, brings about a binding energy of 377.6meV, a value that appears to be within the acceptable range of the adsorption energies for the system under question. Capitalizing on the satisfactory performance of the D3 correction method for similar systems, in a systematic study by Tamijani et al. [79] this approach has been validated and its results have been benchmarked for RG-TiO2(110) utilizing both SAPT(PBE) and experimental adsorption energy data. This scheme, has also led to accurate adsorption energies for methane on rutilelike surfaces as outlined by Weaver [80]. This last study deduces that the D3 dispersioncorrection model along with DFT is sufficiently capable of generating accurate binding energies for physisorbed methane on rutile-like IrO2(110) and RuO2(110) and validates the theoretical binding energetics by their 95% agreement with the experiments. To include the intra-adsorbate interaction effect on the binding strength, adsorption energies at 0.5ML and 1.0ML surface coverages were computed. To do so, 0.5ML of surface coverage was simulated by deposition of CO molecules on every other five-fold Zn atoms of the slab. A monolayer of surface coverage was modelled by placing ad-CO particles onto all of the unsaturated Zn atoms of the surface. Adsorption energies were then calculated. The results are presented in the Table 6. Table 6. Adsorption energies (meV) and, in brackets, equilibrium distances (Å) at 0.5ML and 1.0ML for C-ended isomers.

0.5ML 1.0ML

PBE 236.5[2.53] 149.2[2.57]

optB88-vdW 467.4[2.43] 411.7[2.47]

vdW-DF2 352.3[2.61] 305.4[2.68]

vdW-DF 357.4[2.69] 299.5[2.74]

RPBE-D3 369.4[2.65] 327.7[2.67]

RPBE-D3(BJ) 416.3[2.58] 379.3[2.66]

B3LYP 300 230

LMP2 360 300

Overall, the adsorption energy experiences a shift towards smaller values at higher coverages. Although the distance between any two adjacent adsorbates in 0.5ML of surface coverage is ~6.3Å, the lateral effects are still in-place and lower the adsorption energy by a small amount. As expected, the influence of surface-mediated CO-CO interactions is far more significant for an ad-monolayer. Generally, it is essential to notice that breaking down the interaction energy into adsorbate-adsorbate and adsorbate-adsorbent portions is difficult [6] and can be misleading. The bare lateral interaction energy has been calculated in the past by excluding the surface and solely accounting for the CO-CO interactions using B3LYP. A value of 75meV has been reported for such effects for CO-MgF2(110) system [6]. Based on the assumption that these interactions have a repulsive net impact [81], this value can be added to monolayer adsorption energy to approximate ‘zero-coverage’ limit binding energy. The B3LYP-computed adsorption energy for a monolayer of CO atop MgF2(110) has been previously predicted to be 228meV [6]. As such, the zero-coverage adsorption energy translates to a value as large as 303meV. Compared with highly-accurate binding energy range for CO-MgF2 couple (i.e. 373-385meV), B3LYP is found to underemphasize the adsorption strength. For the same reason, it can be deduced that B3LYP-estimated adsorption energies for 0.5ML and 1.0ML of CO physisorbed onto ZnF2 are slightly underestimated [17]. This is also confirmed by comparing the B3LYP-computed binding energies for CO-ZnF2(110) at 0.5ML and 1.0ML to those of LMP2 one. Clearly, as the surface coverage is increased from half-layer to monolayer, the repulsive intramolecular effects become stronger. Thus, at 1.0 ML of coverage the lateral effects are greater as far as their magnitude. The lateral effects for a halflayer of surface coverage are much smaller and are reported to be between 2-9 meV also using B3LYP [6]. Evidently, increasing the coverage to twice as before (0.5ML→1.0ML) does not double the lateral effects. Moreover, the distance between any two adjacent CO molecules at 0.5ML is longer than 6Å whereas the distance between any two neighboring CO’s at the monolayer level is only slightly greater than 3Å. The ground state intermolecular equilibrium 13

distance for gaseous CO dimer has been measured as 4.0Å and 4.4Å, for C- and O-bonded planar slipped-parallel configurations, respectively [82]. These distances are reported based on broad-band infrared spectroscopic data. As such, monolayer intra-adsorbate repulsive forces are stronger than repulsive interactions at 0.5ML. Our calculations also have indicated the same pattern. Our 0.5ML binding energies are significantly larger than 1.0ML binding energy results. This is true for all the exchange-correlation functionals employed herein. The adsorption energies for 0.5ML of CO adsorbed onto MgF2(110) was previously computed to be 245 meV [6]. This value is also obtained at B3LYP level of theory using valence triple zeta quality Gaussian-type atomic orbitals along with polarization functions as implemented in CRYSTA09 [83,84]. However, the study focusing on this problem also has pointed out that B3LYP underestimates the correlation portion of the adsorption energy by approximately 100 meV [6]. Consequently, a plausible value for the adsorption energy associated with 0.5ML of CO coverage on MgF2 is ~345meV. Given the similarity between MgF2 and ZnF2, the adsorption energy for a half-layer CO-covered ZnF2 must be close to this value. This is also consistent with the LMP2 binding energy calculated by Kaawar and coworkers (i.e. 360meV) [17]. Compared to the reference value, vdW-DF, vdW-DF2 and RPBED3 performance is superlative. The lateral effects associated with the adsorption of CO on surfaces have been formerly studied. As a brief ‘tour d’ horizon’ the combined experimental-theoretical investigation of the lateral interactions between CO adsorbates on Pt(111) can be mentioned. As per Myshlyavtsev and Stishenko account [81], these interactions of CO adsorbates on a metallic substrate are repulsive. Moreover, based on the research conducted by Silvestrelli and co-workers, lateral effects can contribute to a maximum of 25% of the binding energies and they cannot be neglected completely [49]. Our results agree well with both findings. Nowadays, advancements in high-resolution techniques have made possible obtaining real-space images of the molecular layers deposited onto different adsorbents. These experimental methods can be utilized to provide a more realistic picture of the adsorption of CO monolayers atop MgF2 and ZnF2 surfaces. With the experimental data in hand, properties such as binding energies and configurations will be measurable and examining the veracity of the adsorption energies attained via theoretical approaches will become practicable. CO-MgF2(110) interaction energies for 0.5 ML and 1.0 ML of coverage have been formerly calculated by Hammerschmidt and colleagues using LMP2 [5]. Based on their work, at half-coverage limit, the adsorption energy is 314 meV, while with a monolayer of adsorbate, the binding energy becomes 302 meV. Assuming that the CO molecule is to a small degree less penchant for impinging on ZnF2 compared to MgF2, it is evident that optB88-vdW is overpredicting the molecule-substrate interaction strength for both half- and full-coverages. The RPBE-D3(BJ) approach overestimates the binding energy for half-layer and monolayer. The amount of the reduction of the adsorption energy upon increasing the surface coverage has also been measured for CO-TiO2(110) system using temperature-programmed desorption (TPD). Per the study conducted by Dohnálek et al. [85], the CO binding energy at the (110) facet of rutile TiO2 is extrapolated to 435 meV at ‘zero coverage’. This energy decreases to ~348meV, ~198meV and 124meV for 0.5ML, 1.0ML and 1.5ML coverages, respectively. The effect of adsorbate-induced surface relaxations was also examined for the moleculesubstrate couple at 0.1ML, 0.5ML and 1.0ML. In so-doing, both PBE and optB88-vdW methods were employed. Interestingly, it was found that the five-fold Zn atom moves upwards upon accommodating the effect of CO on surface atomic reorganization. This phenomenon has been earlier discovered experimentally and confirmed theoretically for other instances of gassurface interactions. During the CO molecular adsorption at RuO2(110), also a rutile-like substrate, the adsorption site was pulled upwards and in the direction of the adsorbate by 0.2 Å [73]. As another example, dissociative adsorption of CH4 on top of Ir(111) surface can be 14

mentioned, where an Ir metallic atom moved outwards by 0.6 Å [86]. The same has been reported for the sorption of CO atop close-packed transition and noble metal surfaces; an average buckling of the metallic site by 0.1-0.2Å [87]. These instances are interpreted as the outcome of strong attraction between the surface and the adsorbed specimens. This attractive interaction is a result of the superposition of the surface d states at the Fermi energy level and the molecular orbitals of CO, namely, 5σ and 2π* [87]. This interaction draws the surface metallic center towards the adsorbed CO. Very recently, the same structural distortion has been reported in Kaawar et al. for 0.5ML and 1.0ML of CO adsorbates on ZnF2 using B3LYP [17]. The acidic site has been reported to undergo a displacement of 0.03Å. Also, the C≡O distance was reported to be shortened by 0.003-0.006Å. Our results for the amounts by which Zn5f is displaced and C≡O bond is compressed are presented in Tables 7 and 8. Table 7. Zn5f upward displacements upon including the adsorbate-induced relaxation effects (Å).

PBE optB88-vdW B3LYP

0.1ML 0.20 0.15 -

0.5ML 0.17 0.13 -

1.0ML 0.07 0.06 -

0.03

Table 8. C≡O bond length decrease upon including the adsorbate-induced relaxation effects (Å).

PBE optB88-vdW B3LYP

0.1ML 0.008 0.008 -

0.5ML 0.007 0.007 -

1.0ML 0.005 0.006 -

0.003-0.006

Altogether, our estimated atomic displacements for five-fold Zn are overestimated as compared with the B3LYP calculated results. Although, as the surface coverage is increased the agreement between our calculations and reference value is bettered. The PBE-predicted C≡O bond length for free molecule is 1.1355Å. After including the adsorbate-induced relaxation effects, this value decreases by 0.005-0.008Å. The optB88-vdW calculated interatomic distance is 1.1325Å. The bond however is shortened by 0.006-0.008Å as the adsorbate-surface restructuring effects are taken into account. These values are in passable agreement with B3LYP-estimated geometrical predictions. As before, increasing the coverage refines the consistency between our results and literature value. Upon accounting for these effects, the shortest surface-molecule distance is altered too. Our fitting procedure determined an equilibrium distance of 2.52Å for PBE results. After factoring in the impact of adsorbatesurface interplay, the 0.1ML adsorption height was reduced to 2.28 Å. The decline in the molecule-substrate separation is a result of the protrusion of the metallic site from its horizontal plane. The same pattern is observed using the optB88-vdW functional at 0.1ML. The smallest distance between the adsorbate’s C atom and the nearest under-coordinated Zn of the surface changed from 2.57Å to 2.24Å. Including these mutual interactions for 0.1ML increases the optB88-vdW computed adsorption energy to 597.6meV, a clearly overestimated quantity. The PBE-calculated molecule-surface distances before including these effects at 0.5ML and 1.0ML were 2.53Å and 2.57Å, respectively. After including these effects however, the shortest CO-surface distances reduced to 2.31Å and 2.44Å, respectively. Using optB88vdW method, these values were 2.43Å and 2.47Å for 0.5ML and 1.0ML before accounting for adsorbate effects on surface structure. Upon incorporating the adsorbate-surface structural synergic interactions, the minimal energy vertical distances were shortened to 2.26Å and 2.35Å for half- and monolayered surface. Apart from MgF2(110), adsorption of CO onto the (110) plane of other rutile-like surfaces has been examined using both experimental and theoretical approaches. Examples are the sorption of CO atop TiO2(110) [67], SnO2(110) [72] and RuO2(110) [88]. DFT calculations suggested adsorption energies for these systems to be 0.32eV, 0.25eV and 1.2eV, respectively. These energies fall within the scope of previously reported values for dispersion-governed adsorption. DFT predicted metal-C distances for these systems are reported as 2.40 Å, 2.09 Å 15

and 1.95 Å, respectively. The results of the current work communicated that for C-bound isomer and at low coverage, the Zn5f-C distances vary between 2.37Å (calculated via optB88vdW) and 2.69Å (computed via vdW-DF). At higher surface coverages, the shortest moleculeZn5f separation was within 2.43Å (predicted via optB88-vdW at 0.5ML) and 2.74Å (obtained via vdW-DF at 1.0ML). As the surface coverage was increased, the equilibrium distance experienced an elongation as a result of the binding strength weakening. The longest molecule-surface distance at 0.5ML and 1.0ML was predicted through vdW-DF functional. This is commensurate with prior findings that the vdW-DF functional can be too repulsive at short-range [31]. The minimal energy internuclear separation between the C-bridged CO adsorbate and the (110) face of rutile MgF2 is computed to be 2.538 Å using B3LYP [5]. As mentioned earlier, the six-coordinated Mg2+ radius is 0.72 Å. The ionic radius of a six-fold Zn2+ species is 0.74Å. Thus, it is expected for the Zn5f-C bond distance to turn out to be slightly longer than Mg5f-C equilibrium distance. Additionally, the more acidic nature of the Mg centers of MgF2 surface compared to the Zn sites of ZnF2 substrate, dictates that the equilibrium distance for Zn5f-C should be larger than that of Mg5f-C one. As a result, PBE and optB88-vdW both underestimate the equilibrium separation between the adsorbed molecule and the surface.

IV. Concluding remarks and prospective directions To summarize, surface energies and atomic displacements of rutile ZnF2(110) with respect to bulk-truncated structure were predicted using dispersion-free PBE and dispersion-corrected optB88-vdW. The optB88-vdW-relaxed coordinates were then exploited in the calculations. The vdW-DF and vdW-DF2 calculated surface energies were computed at these coordinates and the results were in excellent agreement with literature value of 0.41Jm-2 [3,17]. However, PBE and optB88-vdW under- and overpredict the surface energy. The adsorption of CO onto the (110) face of ZnF2 was simulated using a multitude of dispersion-corrected exchangecorrelation functionals. The vdW-DF and vdW-DF2 computed binding energies were reminiscent of vigorously calculated LMP2 results as reported by Kaawar and colleagues. Although not as accurate as the two aforementioned approaches, zero-damped RPBE-D3 values for 0.5ML and 1.0ML of surface coverage were defensible. Zero-point vibrational frequencies were calculated for free and physisorbed CO using PBE and optB88-vdW methods. Theoretical blue-shifted stretching frequencies were in reasonable accord with experimental data available. The effect of adsorbate-induced surface relaxation was examined and the structural changes were observed in accordance with prior theoretical findings. However, the resultant sorption energy values were overemphasized when including mutual surface-molecule structural interplay. It is difficult to identify a single dispersion-corrected method as the panacea which generates satisfying results for all of the characteristics of interest. Rather, there are groups of vdW-included approaches that can be referred to as reliable and trustworthy for each specific set of properties per system. An in-depth delineation of the small molecule adsorption at rutile-type surfaces is necessary to attain a more realistic picture of the adsorption process. In this respect, newly developed methodologies such as upscaling correction scheme of LMP2→CCSD(T) can be used to reduce the computational effort involved in CCSD(T). This approach has been undertaken for He- and Ar-MgO(100) adsorption systems and has led to plausible adsorption energetics [89,90]. It is safe to say that the modelling procedure of such systems can always be further improved by the accessibility of faster computing processors and the implementation of more computationally efficient codes. The cluster-based results can be utilized to complement periodic DFT-based results. Another major factor in many physical adsorption cases is the effect of competing species. 16

One example of such competition is between any adsorbate and ubiquitous H2O molecules [91]. Given that most of the chemical transformations take place in an aqueous environment, the role of this variable is underscored. The inclusion of the solvation factor will undoubtedly complicate the problem to be attacked. More structurally complicated physisorbates such as pyrrole and CHCl3 (both of which have also been utilized to probe the surface of ZnF2 [1]) can be investigated to foster the physical insight into the adsorption processes on ZnF2 surface. Acknowledgements We wish to extend our sincere gratitude towards Wake Forest University (USA) for allocating extensive computational hardware, as well as the DEAC cluster staffs for their outstanding technical support during this work. Note: The authors of the present work declare no competing financial interest.

References [1] Y. Guo, S. Wuttke, A. Vimont, M. Daturi, J. Lavalley, K. Teinz, E. Kemnitz, J. Mater. Chem. 22 (2012) 14587. [2] I. Suzuki, T. Omata, Y. Shiratsuchi, R. Nakatani, N. Kitamura, S. Otsuka-Yao-Matsuo, Thin Solid Films, 534 (2013) 508. [3] Z. Kaawar, B. Paulus, AIP Conf. Proc., 1653 (2015) 020051. [4] E. Kanaki, G. Sansone, L. Maschio, B. Paulus, Phys. Chem. Chem. Phys. 17 (2015) 18722. [5] L. Hammerschmidt, C. Müller, B. Paulus, J. Chem. Phys. 136 (2012) 124117. [6] Z. Huesges, C. Müller, B. Paulus, C. Hough, N. Harrison, E. Kemnitz, Surf. Sci. 609 (2013) 73. [7] F. Corà, M. Alfredsson, G. Mallia, D.S. Middlemiss, W.C. Mackrodt, R. Dovesi, R. Orlando, Struct. Bond. 113 (2004) 171. [8] M. Darvish Ganji, S.M. Hosseini-Khah, Z. Amini-tabar, Phys. Chem. Chem. Phys. 17 (2015) 2504. [9] F. Chiatti, M. Corno, Y. Sakhno, G. Martra, P. Ugliengo, J. Phys. Chem. C 117 (2013) 25526. [10] L. Díaz Soto, A. Sierraalta, R. Añez, M.A. Chaer Nascimento, J. Phys. Chem. C 119 (2015) 8112. [11] V. Nagarajan, R. Chandiramouli, Appl. Surf. Sci. 385 (2016) 113. [12] R. Blanco, J.M. Orts, Electrochim. Acta 53 (2008) 7796. [13] A. Kiejna, T. Pabisiak, S.W. Gao, J. Phys.: Condens. Matter 18 (2006) 4207. [14] Y. Yamaguchi, K. Tabata, T. Yashima, J. Mol. Struct. (Theochem) 714 (2005) 221. [15] H.H. Klauss, N. Sahoo, P.C. Kelires, T.P. Das, W. Potzel, M. Kalvius, M. Frank, W. Kreische, Hyperfine Interact. 60 (1990) 853. 17

[16] Y. Wan-Lun, Z. Min-Guang, J. Phys. C: Solid State Phys. 18 (1985) L1087. [17] Z. Kaawar, C. Müller, B. Paulus, Surf. Sci. 656 (2017) 48. [18] P. Zeppenfeld, Surface and Interface Science: Solid-Gas Interfaces, first ed., Edited by: K. Wandelt, Wiley-VCH Verlag GmbH & Co. KGaA, 2015. [19] T. Trevethan, A.L. Shluger, J. Phys. Chem. C 111 (2007) 15375. [20] A.D. Boese, J. Sauer, J. Comput. Chem. 37 (2016) 2374. [21] C.A. Downing, A.A. Sokol, C.R.A. Catlow, Phys. Chem. Chem. Phys. 16 (2014) 184. [22] T. Krahl, A. Vimont, G. Eltanany, M. Daturi, E. Kemnitz, J. Phys. Chem. 111 (2007) 18317. [23] M. Buchholz, X. Yu, C. Yang, S. Heißler, A. Nafedov, Y. Wang, C. Wöll, Surf. Sci. 652 (2016) 247. [24] A.D. Becke, J. Chem. Phys. 140 (2014) 18A301. [25] A.E. Mattsson, P.A. Schultz, M.P. Desjarlais, T.R. Mattsson, K. Leung, Modelling Simul. Mater. Sci. Eng. 13 (2005) R1. [26] S. Grimme, S. Ehlrich, L. Goerigk, J. Comput. Chem. 32 (2011) 1456. [27] E.R. Johnson, A.D. Becke, J. Chem. Phys. 124 (2006) 174104. [28] C. Adamo, V. Barone, J. Chem. Phys. 110 (1999) 6158. [29] J. Klimeš, A. Michaelides, J. Chem. Phys. 137 (2012) 120901. [30] M. Dion, H. Rydberg, E. Schröder, D.C. Langreth, B.I. Lundqvist, Phys. Rev. Lett. 92 (2004) 246401. [31] J. Klimeš, D.R. Bowler, A. Michaelides, Phys. Rev. B. 83 (2011) 195131. [32] K. Yang, J. Zheng, Y. Zhao, D.G. Truhlar, J. Chem. Phys. 132 (2010) 164117. [33] K. Lee, È.D. Murray, L. Kong, B.I. Lundqvist, D.C. Langreth, Phys. Rev. B. 82 (2010) 081101. [34] K. Lee, Y. Morikawa, D.C. Langreth, Phys. Rev. B. 82 (2010) 155461. [35] X. Ma, A. Genest, L. Spanu, N. Rösch, Comput. Theor. Chem. 1069 (2015) 147. [36] K. Berland, V.R. Cooper, K. Lee, E. Schröder, T. Thonhauser, P. Hyldgaard, B.I. Lundqvist, Rep. Prog. Phys. 78 (2015) 066501. [37] J. Carrasco, W. Liu, A. Michaelides, A. Tkatchenko, J. Chem. Phys. 140 (2014) 084704. [38] M.J. Tillotson, P.M. Brett, R.A. Bennett, R. Grau-Crespo, Surf. Sci. 632 (2015) 142. [39] G. Kresse, J. Furthmüller, Phys. Rev. B 54 (1996) 11169. [40] G. Kresse, J. Furthmüller, Comput. Mater. Sci. 6 (1996) 15. 18

[41] H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13 (1976) 5188. [42] J. Haines, J. Léger, F. Gorelli, D. Klug, J.S. Tse, Z. Li, Phys. Rev. B 64 (2001) 134110. [43] S.H. Lin, J.L. Kuo, Phys. Chem. Chem. Phys. 16 (2014) 20763. [44] J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [45] P. Zhao, Y. He, D. Cao, X. Wen, H. Xiang, Y. Li, J. Wang, H. Jiao, Phys. Chem. Chem. Phys. 17 (2015) 19446. [46] M. Bajdich, J.K. Nørskov, A. Vojvodic, Phys. Rev. B 91 (2015) 155401. [47] G. Makov, M.C. Payne, Phys. Rev. B 51 (1994) 4014. [48] J. Tao, A.M. Rappe, Phys. Rev. Lett. 112 (2014) 106101. [49] P.L. Silvestrelli, A. Ambrosetti, 185 (2016) 1. [50] W.M. Haynes, Editor-in-chief. CRC Handbook of Chemistry and Physics. CRC PRESS, 86th edition, 2015-2016. [51] National Institute of Standard and Technology (2011) NIST Chemistry WebBook. http://webbook.nist.gov [52] F. Tran, R. Laskowski, P. Blaha, K. Schwarz, Phys. Rev. B 75 (2007) 115131. [53] K. Chuasiripattana, O. Waarschkov, B. Delley, C. Stampfl, Surf. Sci. 604 (2010) 1742. [54] W.G. Reimers, M.A. Baltanás, M.M. Branda, J. Mol. Model. 20 (2014) 2270. [55] R. Dawes, X.G. Wang, T. Carrington Jr., J. Phys. Chem. A 117 (2013) 7612. [56] T.R. Esch, I. Gadaczek, T. Bredow, Appl. Surf. Sci. 288 (2014) 275. [57] T.A. Mellan, K.P. Maenetja, P.E. Ngoepe, S.M. Woodley, C. Richard, A. Catlow, R. Grau-Crespo, J. Mater. Chem A. 1 (2013) 14879. [58] L. Hammerschmidt, L. Maschio, C. Müller, B. Paulus, J. Chem. Theory Comput. 11 (2015) 252. [59] W. Wei, Y. Dai, B. Huang, J. Phys. Chem. C 115 (2011) 18597. [60] R.D. Shannon, Acta Crystallogr. A32 (1976) 751. [61] Y.Q. Jia, J. Solid State Chem. 95 (1991) 184. [62] T.M. Muscenti, G.V. Gibbs, D.F. Cox, Surf. Sci. 594 (2005) 70. [63] D.B. Rogers, R.D. Shannon, A.W. Sleight, J.L. Gillson, Inorg. Chem. 8 (1969) 841. [64] D. Stodt, H. Noei, C. Hätting, Y. Wang, Phys. Chem. Chem. Phys. 15 (2013) 466. [65] F. Munoz, A.H. Romero, J. Mejía-López, I.V. Roshchin, R.I. González, M. Kiwi, J. Magn. Magn. Mater. 393 (2015) 226. 19

[66] D.J. Mowbray, J.I. Martínez, F. Calle-Vallejo, J. Rossmeisl, K.S. Thygesen, K.W. Jacobsen, J.K. Nørskov, J. Phys. Chem. C 115 (2011) 2244. [67] J. Blomqvist, L. Lehman, P. Salo, Phys. Status Solidi B 249 (2012) 1046. [68] T. Diemant, H. Hartmann, J. Bansmann, R.J. Behm, Surf. Sci. 652 (2016) 123. [69] M. Pykavy, V. Staemmler, O. Seiterth, H.-J. Freund, Surf. Sci. 479 (2001) 11. [70] N. Vulliermet, T.A. Wesolowski, J. Weber, Collect. Czech. Chem. Commun. 63 (1998) 1447. [71] N.G. Petrik, G.A. Kimmel, Phys. Chem. Chem. Phys. 16 (2014) 2338. [72] M. Melle-Franco, G. Pacchioni, Surf. Sci. 461 (2000) 54. [73] A.P. Seitsonen, Y.D. Kim, M. Knapp, S. Wendt, H. Over, Phys. Rev. B 65 (2001) 035413. [74] C.L. Angell, P.C. Schaffer, J. Phys. Chem. 70 (1966) 1413. [75] W.M. Huo, J. Chem. Phys. 43 (1965) 624. [76] M. Setvin, M. Buchholz, W. Hou, C. Zhang, B. Stöger, J. Hulva, T. Simschitz, X. Shi, J. Pavelec, G.S. Parkinson, M. Xu, Y. Wang, M. Schmid, C. Wöll, A. Selloni, U. Diebold, J. Phys. Chem. C 119 (2015) 21044. [77] C. Lemire, R. Meyer, V.E. Henrich, Sh. Shaikhutdinov, H.-J. Freund, 572 (2004) 103. [78] Í.P. de Lima, J.R. dos S. Politi, R. Gargano, J.B.L. Martins, Theor. Chem. Acc. 134 (2015) 1. [79] A.A. Tamijani, A. Salam, M.P. de Lara-Castells, J. Phys. Chem. C 120 (2016) 18126. [80] J.F. Weaver, Chem. Rev. 113 (2013) 4164. [81] A.V. Myshlyavtsev, P.V. Stishenko, Surf. Sci. 642 (2015) 51. [82] M. Rezaei, S. Sheybani-Deloui, N. Moazzen-Ahmadi, K.H. Michaelian, A.R.W. McKellar, J. Phys. Chem. A 117 (2013) 9612. [83] R. Dovesi, R. Orlando, B. Civalleri, R. Roetti, V.R. Saunders, C.M. Zicovich-Wilson, Z. Kristallogr. 220 (2005) 571. [84] R. Dovesi, V.R. Saunders, R. Roetti, R. Orlando, C.M. Zicovich-Wilson, F. Pascale, B. Civalleri, K. Doll, N.M. Harrison, I.J. Bush, P.D’Arco, M. Llunell, CRYSTAL09 User’s Manual, University of Torino, 2009. [85] Z. Dohnálek, J. Kim, O. Bondarchuk, J.M. White, B.D. Kay, J. Phys. Chem. B 110 (2006) 6229. [86] G. Henkelman, H. Jónsson, Phys. Rev. Lett. 86 (2001) 664. [87] M. Gajdoš, A. Eichler, J. Hafner, J. Phys.: Condens. Matter 16 (2004) 1141. [88] Y.D. Kim, A.P. Seitsonen, H. Over, Phys. Rev. B 63 (2001) 115419. 20

[89] R. Martinez-Casado, D. Usvyat, L. Maschio, G. Mallia, S. Casassa, J. Ellis, M. Schütz, N.M. Harrison, Phys. Rev. B 89 (2014) 205138. [90] D. Usvyat, K. Sadeghian, L. Maschio, M. Schütz, Phys. Rev. B 86 (2012) 045412. [91] M. Patel, F.F. Sanches, G. Mallia, N.M. Harrison, Phys. Chem. Chem. Phys. 16 (2014) 21002.

21

Two-fold F

Three-fold F Five-fold Zn

Six-fold Zn

Figure 1. The optB88-vdW-relaxed monolayer of CO adsorbed atop rutile ZnF2(110) surface, as viewed along [100] (top) and [010] (bottom) crystallographic axes. Turquoise, purple, black and red spheres are Zn, F, C and O, respectively.

22

Highlights Atomic displacements of the (110) plane of rutile-structured ZnF2 are reported. The vdW-DF and vdW-DF2 computed surface energies of r-ZnF2 afford excellent agreement with that of B3LYP one. Binding energies of CO at (110) facet of r-ZnF2 at 0.1ML, 0.5ML and 1.0ML of surface coverage are calculated using a manifold of dispersion-corrected DFT methods. At 0.5ML and 1.0ML of coverage limits, the vdW-DF and vdW-DF2 produce adsorption energies as accurately as the less computationally efficient LMP2 and more precisely than dispersion-less B3LYP.

graphical abstract

23