CO2 stability on the Ni low-index surfaces: van der Waals corrected DFT analysis

Catalysis Communications 80 (2016) 33–38

Contents lists available at ScienceDirect

Catalysis Communications journal homepage: www.elsevier.com/locate/catcom

Short communication

CO2 stability on the Ni low-index surfaces: van der Waals corrected DFT analysis Kamil Czelej ⁎, Karol Cwieka, Krzysztof Jan Kurzydlowski Faculty of Materials Science and Engineering, Warsaw University of Technology, Woloska 141, 02-507 Warsaw, Poland

a r t i c l e

i n f o

Article history: Received 20 January 2016 Received in revised form 6 March 2016 Accepted 24 March 2016 Available online 31 March 2016 Keywords: Heterogeneous catalysis Adsorption Decomposition Nickel CO2 DFT

a b s t r a c t Carbon dioxide stability on the nickel low-index surfaces has been studied by means of van der Waals corrected spin-polarized density functional theory. A number of possible CO2/Nisurface conformations with negative adsorption energy were identified. The partial density of states combined with the effective bond order results indicate significant activation of the C_O bond by enhanced charge transfer and shift of the antibonding molecular orbital below the Fermi level. On the basis of the potential energy diagrams, high mobility of COδ− 2 moiety on the Ni low-index surfaces and thermodynamic preference for decomposition of CO2 to surface bound CO and O were predicted. The Ni(100) surface was found to be the most efficient in terms of CO2 conversion to CO and O. © 2016 Elsevier B.V. All rights reserved.

1. Introduction The greenhouse effect, due to rapid increase of carbon dioxide (CO2) emission into the atmosphere, has become one of the major concerns of modern societies. It is estimated that CO2 contribution to the global warming has exceeded 60% [1] and its atmospheric concentration keeps rising at an accelerating rate [2]. For this reason, significant effort has been devoted to investigating efficient ways to reduce CO2 and convert it into useful organic molecules, as a precursor to fuels and chemicals feedstocks [3–8]. In particular, the heterogeneous catalytic conversion of CO2 on noble metals and transition metal based alloys has attracted immense scientific interest. One of the most relevant problems associated with heterogeneous catalytic conversion of CO2 is the formation of coke on the surface, which in turn leads to deactivation of the catalysts and limits their lifetime. Although noble metal based catalysts are effective and at the same time less sensitive to coking than transition metal based catalysts, their availability is limited and cost is very high. Therefore, the development of transition metal based catalysts of high efficiency, resistant to coking seems economically justified. Amongst transition metals nickel is cheap, available and has been successfully applied in different fields, such as catalytic conversion of CO2 to methane [9,10], catalytic conversion of CO2 to CO with CH4 [11] or base material for electrodes in molten carbonate fuel cell (MCFC) [12]. In order to improve efficiency of Ni based materials towards CO2 conversion and design a new Ni based catalysts, an elemental ⁎ Corresponding author. E-mail address: [email protected] (K. Czelej).

http://dx.doi.org/10.1016/j.catcom.2016.03.017 1566-7367/© 2016 Elsevier B.V. All rights reserved.

catalytic steps, such as CO2 adsorption and decomposition on Ni surfaces, have to be clarified. Over the past three decades, the development of computational ab initio methods allowed to investigate an interaction between CO2 molecule and the Ni low-index surfaces [13–20]. Both chemisorbed and physisorbed species were predicted via DFT results. The electronic structure calculation of the chemisorbed conformations revealed the activation of the C_O bond and charge transfer from the Ni surface to the molecule, giving rise to formation of a COδ− moiety [13–15]. 2 Wang et al. [14] carried out the bond nature analysis of the adsorbed CO2 and concluded that chemisorption ability decreased in the order Ni(110) N Ni(100) N Ni(111). Recently, Ding et al. [17] and Dri et al. [18] confirmed the coexistence of a number of chemisorbed CO2 conformations on Ni(110) using a combination of DFT calculations with experimental techniques, such as scanning tunneling microscopy, photoelectron spectroscopy and high-resolution electron energy loss spectroscopy. Decomposition of CO2 to CO and O, has been investigated on Ni(100) and Ni(110). Liu et al. [19] determined the minimum energy pathway for CO2 decomposition on a Ni(100) surface and calculated the total energy barrier of 11.1 kcal/mol (~0.48 eV). In our recent study [20] we carried out a set of nudged elastic band (NEB) calculations for CO2 decomposition on Ni(110) in order to clarify the reaction mechanism. It was concluded that breakage of the coordinated C\\O bond is preceded by surface diffusion to the one energetically favored conformation and the total energy barrier was found to be about 0.44 eV. In the present study, the stability of CO2 on the Ni low-index surfaces was investigated in the framework of spin-polarized DFT calculations. An impact of van der Waals interaction on the adsorption energy and

34

K. Czelej et al. / Catalysis Communications 80 (2016) 33–38

Fig. 1. Top view of possible adsorption sites on Ni(100), Ni(110) and Ni(111).

the geometry of CO2 species was shown. The transition state theory calculations within the nudged elastic band and dimer formalism were further utilized to find relevant transition states and calculate diffusion and decomposition energy barriers. 2. Computational method The structural relaxation and total energy calculations were performed in the framework of spin-polarized density functional theory (SP-DFT) as implemented in the Vienna Ab initio Simulation Package (VASP) [21]. The Perdew–Burke–Ernzerhof (PBE) general gradient approximation (GGA) functional [22] was used to estimate the exchange-correlation energy. The Kohn–Sham one-electron states were expanded in a plane wave basis set, whereas ionic cores were approximated by the projector augmented wave (PAW) method [23]. The description of nonlocal long-range interactions between adsorbed molecules and metal surfaces was improved by employment of the vdW correction proposed by S. Grimme (DFT-D2) [24]. Dipole correction was taken into account for all calculations [25]. In order to simulate an individual CO2 molecule on Ni(100), Ni(110) and Ni(111) surfaces, a p(3 × 3) slab geometries, corresponding to a 1/9 ML coverage, containing 5 metal layers were constructed from the optimized bulk unit cell of nickel. The calculated equilibrium lattice constant of 3.514 Å and magnetic moment of 0.63 μB are in good agreement with the experimental values of 3.516 Å and 0.61 μB, respectively [26]. The p(3 × 3) slab supercell with 5 nickel layers and a vacuum level of ~ 15 Å was large enough to avoid the spurious interaction between CO2 molecule and its periodic images. Test calculations for a larger p(4 × 4) supercell gave a numerical difference in CO2 adsorption energy of less than 0.01 eV. Convergence of the total energy for the slabmolecule system was assured by setting a 450 eV cutoff energy, and the Brillouin zone sampling in a 5 × 5 × 1, 4 × 5 × 1 and 5 × 5 × 1 Monkhorst–Pack k-point set [27] for Ni(100), Ni(110) and Ni(111), respectively. For geometry optimization, the top 3 layers of nickel atoms together with the CO2 molecule were allowed to relax towards a minimum of the total energy, while the bottom 2 layers were kept frozen. After all forces acting on atoms dropped below 0.01 eV/Å and energy convergence criterion of 1e − 5 eV was reached, the relaxations were terminated. A number of inequivalent adsorption sites for the Ni low-index surfaces can be distinguished as shown in Fig. 1. To determine all stable conformations, plenty of different initial configurations, corresponding to various molecular orientations for each adsorption site were examined. The adsorption energy per CO2 molecule is given by: ΔEads ¼ EðCO2 =Nislab Þ−EðNislab Þ−EðCO2 Þ

where the first term is the computed total energy for the slab with adsorbed CO2 on the Ni surface, the second term is the energy of the bare metal surface and the third term is the total energy of free CO2. Therefore, exothermic adsorption corresponds to a negative Δ Eads value while endothermic adsorption corresponds to a positive Δ Eads value. Base on Δ Eads, the relative stability of different conformations can be analyzed. In order to quantify charge transfer from the Ni low-index surfaces to the molecule, a Bader analysis was carried out for all CO2/Nisurface conformations using the code developed by Henkelman et al. [28]. Effective bond order (EBO) calculations combined with partial density of states (PDOS) results were utilized to investigate the bonding nature of CO2 molecule on the Ni(100), Ni(110) and Ni(111) surfaces. Transition state searches were performed using a combination of climbing image nudged elastic band (CI-NEB) method [29] and dimer method. The NEB method allows one to determine a reaction pathway when the initial and final configurations are known. By contrast, the dimer method requires only one initial configuration and the initial direction along the dimer, which can be determined, for instance, from pre-converged NEB pathway. To achieve force (0.03 eV/Å) and energy (1e− 5 eV) convergence criteria, global force NEB optimizations were

Fig. 2. Possible adsorption conformations of CO2 on the nickel low-index surfaces.

K. Czelej et al. / Catalysis Communications 80 (2016) 33–38

35

Table 1 Computed geometrical parameters and chemisorption energies of CO2 species on the Ni low-index surfaces. Surface

100

110a

111

a

Conformation

H-4f H-3f H-2f B-2f H-3f SB-4f H-5f LB-2f SB-2f Phys. Phys. B-2f H-4f H-3f B-4f

Structure geometry dO(1)\ \C, Å

dO(2)\ \C, Å

dO(1)\ \Ni, Å

dO(2)\ \Ni, Å

dC\ \Ni, Å

α(O\ \C\ \O), deg

1.383 1.274 1.311 1.277 1.284 1.246 1.306 1.258 1.370 1.178 1.178 1.278 1.280 1.341 1.241

1.226 1.274 1.228 1.219 1.286 1.246 1.306 1.240 1.242 1.178 1.178 1.216 1.278 1.206 1.241

1.958 2.001 1.937 1.947 1.947 2.033 2.064 1.972 2.002 2.836 3.377 1.967 2.094 2.044 2.100

2.786 2.001 2.481 2.642 1.971 2.040 2.065 2.222 2.563 2.838 3.377 2.654 2.097 2.740 2.100

2.020 1.897 1.914 1.925 1.916 2.018 1.898 1.889 2.037 3.714 3.136 1.956 1.993 1.903 2.096

123.5 128.7 129.2 135.6 124.4 138.9 121.6 136.3 122.1 178.7 179.4 136.0 131.9 129.6 140.1

ΔEads, eV

Bader charge of CO2 moiety

−0.65 −0.50 −0.39 −0.38 −0.77 −0.68 −0.67 −0.66 −0.55 −0.30 −0.22 −0.18 −0.16 −0.17 −0.11

−0.849 −0.730 −0.728 −0.546 −0.832 −0.606 −0.978 −0.593 −0.935 0.000 0.000 −0.518 −0.742 −0.618 −0.569

Values taken from our previous work [20].

performed using fast inertial relaxation engine and limited-memory Broyden–Fletcher–Goldfarb–Shanno ionic optimizers [30,31]. The transition states were further verified by vibration analysis, carried out within the dynamical matrix code. 3. Results and discussion 3.1. CO2 adsorption on the Ni low-index surfaces Interaction of CO2 with the Ni low-index surfaces gives rise to formation of both physisorbed and chemisorbed conformations. In case of physisorption, the molecule is weakly bound to the Ni surfaces via van der Waals forces with nearly unchanged linear geometry. Only two physisorbed species on Ni(110) and Ni(111) surfaces were found (see

Fig. 2). In general, DFT does not properly account for the van der Waals interactions, responsible for binding molecular physisorbed states on surfaces. Therefore, a semi-empirical addition of dispersive forces to GGA functionals might be crucial to suitably describe the geometry and the adsorption energy of physisorbed species. In fact, the DFT-D2 results indicate a significant reduction in the distance between CO2 molecule and Ni surface by about 0.73 Å for Ni(110) and 0.65 Å for Ni(111) as well as higher adsorption energies by about 0.41 and 0.42 eV for Ni(110) and Ni(111), respectively, in comparison to GGA results (see Supplementary Material). In case of chemisorption, however, the chemical bond between CO2 molecule and Ni surface atoms is created, causing changes in electronic structure and geometry of the molecule. A number of chemisorbed CO2 conformations on Ni(100), Ni(110) and Ni(111) were found, as shown in Fig. 2 and Supplementary Material (Fig. S1). The computed structural parameters of the chemisorbed species together with adsorption energies are summarized in Table 1. An influence of van der Waals interaction on the structural parameters of the chemisorbed species is negligibly small, however, the adsorption energy changes significantly by about 0.39, 0.41 and 0.42 eV for Ni(100), Ni(110) and Ni(111), respectively (see Supplementary Material). It is noteworthy that for the Ni(111) surface the van der Waals contribution determines stability of CO2 adsorption by changing its character from endothermic to exothermic. Amongst 6 stable CO2 conformations found on Ni(100) the most stable one is H-4f where CO2 occupies a hollow site interacting with four top-layer Ni atoms via one C_O bond. For Ni(110), five stable chemisorbed CO2 Table 2 Computed effective bond orders for CO2/Nisurface conformations. Surface

100

110a

111

Fig. 3. Partial density of states (PDOS) projected onto molecular orbitals for isolated CO2 molecule and CO2 chemisorbed on the Ni low-index surfaces.

a

Conformation

Effective bond order O(1)\ \C

O(2)\ \C

O(1)\ \Ni

O(2)\ \Ni

C\ \Ni

CO2

2.00

2.00

–

–

–

H-4f H-5f H-2f B = 2f H-3f SB-4f SB-2f H-5f LB-2f H-2f B-2f H-4f B-4f H-3f

1.34 1.73 1.60 1.74 1.69 1.85 1.77 1.53 1.82 1.40 1.73 1.69 1.90 1.44

2.13 1.73 2.09 2.14 1.68 1.84 2.00 1.53 1.94 2.00 2.15 1.71 1.91 2.21

0.52 0.48 0.58 0.57 0.57 0.42 0.54 0.38 0.55 0.47 0.53 0.39 0.35 0.41

0.07 0.48 0.13 0.08 0.52 0.40 0.13 0.38 0.24 0.11 0.08 0.38 0.35 0.05

0.50 0.73 0.66 0.68 0.73 0.43 0.70 0.70 0.69 0.46 0.62 0.50 0.34 0.75

Values taken from our previous work [20].

36

K. Czelej et al. / Catalysis Communications 80 (2016) 33–38

conformations have been found with the most stable one assigned as a H-3f, where CO2 molecule is bound to three top-layer Ni atoms in an asymmetric way and the molecular plane is tilted at an angle of 28° from the metal surface. For Ni(111), four stable CO2 conformations with very similar adsorption energies have been revealed. The most stable one is B-2f, where CO2 molecule is placed in a bridge site. On the basis of the computed adsorption energies, the ability of CO2

chemisorption on the Ni low-index surfaces in the order Ni(110) N Ni(100) N Ni(111) has been confirmed. According to the van der Waals corrected DFT results, CO2 chemisorption on the Ni low-index surfaces is exothermic and differences between the chemisorption energies of the investigated species are small. Therefore, one can expect a coexistence of different CO2 conformations at elevated temperatures.

Fig. 4. Potential Energy diagrams for different CO2 adsorbed conformations on the Ni low-index surfaces with relevant diffusion and decomposition barriers: a) Ni(100), b) Ni(110), c) Ni(111). The zero energy is set to the isolated CO2 molecule at infinite distance from the surface.

K. Czelej et al. / Catalysis Communications 80 (2016) 33–38

37

In order to clarify the bonding nature of CO2/Nisurface system, the partial density of states (PDOS) projected onto molecular orbitals together with the effective bond order calculation results were investigated. The PDOS for all chemisorbed species are very similar therefore, only a selection is displayed on Fig. 3. In isolated CO2 the highest occupied molecular orbital (HOMO) is the degenerate 1πg bonding orbital, while the lowest unoccupied molecular orbital (LUMO) is the degenerate 2πu antibonding orbital. When chemisorption takes place, both, 1πg and 2πu bands drop below the Fermi level and hybridize with d bands of Ni. As a result, a chemical bond between the molecule and the Ni surface is created. The effective bond order calculation results (Table 2) indicate that Ni\\C bond is mainly responsible for the binding of the CO2 molecule to the Ni low-index surfaces. Due to the charge transfer from the Ni d orbitals into the molecular 2πu antibonding orbital (see the Bader analysis results in Table 1), the CO2 molecule bends, forming a negatively charged moiety on the Ni surface. The stronger the charge transfer, the higher deterioration of the C\\O bond leading to its activation.

even changes its character from endothermic to exothermic. Therefore, dispersion interaction cannot be neglected in anticipation of stability of adsorbed species. Taking into account relatively low diffusion energy barriers, high mobility of COδ− 2 moiety can be expected. A negative driving force of the CO2/Nisurface → CO/Nisurface + O/Nisurface reaction makes the decomposition process thermodynamically favorable. On the basis of the computed reaction energies and decomposition energy barriers, the Ni(100) enables the most efficient CO2 conversion yield amongst Ni low-index surfaces.

3.2. Surface diffusion and decomposition of CO2 on the Ni low-index surfaces

Appendix A. Supplementary data

Due to the thermal activation of chemisorbed conformations, the Nisurface/CO2 system may transform, undergoing various transition states. The stability of particular CO2 conformation on the nickel surface depends on a number of factors, such as molecule–metal interaction strength, mobility of chemisorbed species, activation degree of the C\\O bond, etc. In order to shed some light on possible transformations within the Nisurface/CO2 system, the potential energy diagrams for different CO2 adsorbed conformations on the Ni low-index surfaces have been constructed, as shown in Fig. 4. All transition states were verified by vibration analysis; as expected, each transition state is characterized by one imaginary mode, confirming the calculated transitions states (see Supplementary Material). Additionally, detailed geometries of the transition states were shown in Supplementary Material. As can be seen in the energy diagrams, decomposition of CO2 to CO and O is thermodynamically favored due to the negative value of the CO2/ Nisurface → CO/Nisurface + O/Nisurface reaction energy. The driving force of this reaction measured as the difference between the total energy of the (C + O)/Nisurface and CO2/Nisurface conformations decreases in the order Ni(100) N Ni(111) N Ni(110). A similar conclusion was drawn by Wang et al. [13]. The reaction rate however, is a function of total energy barrier in given temperature. According to the results presented in Fig. 4, the CO2 decomposition reaction rate on the Ni low-index surfaces is in the order of Ni(100) N Ni(111) ~ Ni(110). As shown in previous paragraph, the chemisorption energies of the investigated species do not vary significantly for a given Ni surface and a number of chemisorbed species may coexist at elevated temperatures. The calculated surface diffusion barriers are relatively low and therefore, high mobility of COδ− 2 moiety can be expected. For all investigated Ni low-index surfaces the decomposition reaction mechanism is quite similar. The breaking of the coordinated C\\O bond is preceded by the surface diffusion of CO2 moiety to the chair-like geometry in the hollow site. Detailed explanation of this phenomenon for CO2/Ni(110) system can be found in our previous work [20]. 4. Conclusions CO2 stability on the Ni low-index surfaces was examined using the van der Waals corrected SP-DFT method. A number of inequivalent chemisorption sites, which might be simultaneously populated by CO2 at thermal equilibrium, have been revealed. The relative stability of all investigated conformations together with the ability of CO2 chemisorption on the Ni low-index surfaces in the order Ni(110) N Ni(100) N Ni(111) confirm previous theoretical and experimental observations. It has been found that the van der Waals contribution to the total energy significantly increases the adsorption energy, and for the Ni(111) surface

Acknowledgments Computing resources were provided by High Performance Computing facilities of the Interdisciplinary Centre for Mathematical and Computational Modeling (ICM) of the University of Warsaw under Grant No. G61-4.

Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.catcom.2016.03.017. References [1] A. Albo, P. Luis, A. Irabin, Carbon dioxide capture from flue gases using a cross-flow membrane, 49 (2010) 11045–11051. [2] S. Kaneco, M. Yabuuchi, H. Katsumata, T. Suzuki, K. Ohta, ACS Division of Fuel, Chemistry 47 (2002) 71. [3] A.J. Morris, G.J. Meyer, E. Fujita, Molecular approaches to the photocatalytic reduction of carbon dioxide for solar fuels, Acc. Chem. Res. 42 (12) (2009) 1983–1984. [4] G. Centi, S. Perathoner, Catalysis: role and challenges for a sustainable energy, Top. Catal. 52 (8) (2009) 948–961. [5] S.C. Roy, O.K. Varghese, M. Paulose, Toward solar fuels: photocatalytic conversion of carbon dioxide to hydrocarbons, ACS Nano 4 (3) (2010) 1259–1278. [6] J. Baltrusaitis, V.H. Grassian, Surface reactions of carbon dioxide at the adsorbed water–iron oxide interface, J. Phys. Chem. B 109 (25) (2005) 12227–12230. [7] R.D. Richardson, E.J. Holland, B.K. Carpenter, A renewable amine for photochemical reduction of CO2, Nat. Chem. 3 (4) (2011) 301–303. [8] V. Havran, M.P. Duduković, C.S. Lo, Conversion of methane and carbon dioxide to higher value products, Ind. Eng. Chem. Res. 50 (12) (2011) 7089–7100. [9] A. Marinoiu, C. Cobzaru, M. Raceanu, M. Varlam, E. Carcadea, C. Cernatescu, I. Stefanescu, Carbon dioxide conversion to methane over supported nickel base catalysts, Rev. Roum. Chim. 60 (2) (2015) 249–256. [10] Y. Liu, Z. He, L. Zhou, Z. Hou, W. Eli, Simultaneous oxidative conversion and CO2 reforming of methane to syngas over Ni/vermiculite catalysts, Catal. Commun. 42 (2013) 40–44. [11] H. Koh, S. Lee, J. Choi, Conversion of CO2 to CO with CH4 over Ni/SiO2 catalyst: photoacustic measurements of reaction rate, Bull. Kor. Chem. Soc. 25 (8) (2004) 1253–1256. [12] E. Antolini, The stability of molten carbonate fuel cell electrodes: a review of recent improvements, Appl. Energy 88 (2011) 4274–4293. [13] S.G. Wang, D.B. Cao, Y.W. Li, J. Wang, H. Jiao, Chemisorption of CO2 on nickel surfaces, J. Phys. Chem. B 109 (40) (2005) 18956–18963. [14] S.G. Wang, X.Y. Liao, D.B. Cao, C.F. Huo, Y.W. Li, J. Wang, H. Jiao, Factors controlling the interaction of CO2 with transition metal surfaces, J. Phys. Chem. C 111 (2007) 16934–16940. [15] H.J. Freund, R.P. Messmer, On the bonding and reactivity of CO2 on metal surfaces, Surf. Sci. 172 (1) (1986) 1–30. [16] X. Ding, L. De Rogatis, E. Vesselli, A. Baraldi, G. Comelli, R. Rosei, L. Savio, L. Vattuone, M. Rocca, P. Fornasiero, F. Ancilotto, A. Baldereschi, M. Peressi, Interaction of carbon dioxide with Ni(110): a combined experimental and theoretical study, Phys. Rev. B 76 (2007) 195425. [17] X. Ding, V. Pagan, M. Peressi, F. Ancilotto, Modeling adsorption of CO2 on Ni(110) surface, Mater. Sci. Eng. C 27 (5–8) (2007) 1355–1359. [18] C. Dri, A. Peronio, E. Vesselli, C. Africh, M. Rizzi, A. Baldereschi, M. Peressi, G. Comelli, Imaging and characterization of activated CO2 species on Ni(110), Phys. Rev. B 82 (2010) 165403. [19] C. Liu, T.R. Cundari, A.K. Wilson, CO2 reduction on transition metal (Fe, Co, Ni, and Cu) surfaces: in comparison with homogeneous catalysis, J. Phys. Chem. C 116 (February 2012) 5681–5688. [20] K. Czelej, K. Cwieka, T. Wejrzanowski, P. Spiewak, K.J. Kurzydlowski, Decomposition of activated CO2 species on, Catal. Commun. 74 (2016) 65–70. [21] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (October 1996) 11169–11186. [22] J.P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77 (1996) 3865–3868. [23] P.E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50 (December 1994) 17953–17979.

38

K. Czelej et al. / Catalysis Communications 80 (2016) 33–38

[24] S. Grimme, Semiempirical GGA-type density functional constructed with a longrange dispersion correction, J. Comput. Chem. 27 (15) (2006) 1787–1799. [25] L. Bengtsson, Dipole correction for surface supercell calculations, Phys. Rev. B 59 (19) (1999) 12301–12304. [26] R.T. Downs, M. Hall-Wallace, The American Mineralogist crystal structure database, Am. Mineral. 88 (2003) 247–250. [27] H.J. Monkhorst, J.D. Pack, Special points for Brillouin zone integrations, Phys. Rev. B 13 (12) (January 1976) 5188–5192. [28] G. Henkelman, A. Arnaldsson, H. Jonsson, A fast and robust algorithm for Bader decomposition of charge density, Comput. Mater. Sci. 36 (2006) 354–360.

[29] G. Henkelmann, B.P. Uberuaga, H. Jónsson, A climbing image nudged elastic band method for finding saddle points and minimum energy paths, J. Chem. Phys. 113 (22) (2000) 9901–9904. [30] D. Sheppard, R. Terrell, G. Henkelman, Optimization methods for finding minimum energy paths, J. Chem. Phys. 128 (2008) 134106. [31] D. Sheppard, G. Henkelman, Paths to which the nudged elastic band converges, J. Comput. Chem. 32 (2011) 1769–1771.