Comprehensive Phase Diagrams of MoS2 Edge Sites Using

Subscriber access provided by TUFTS UNIV

C: Surfaces, Interfaces, Porous Materials, and Catalysis 2

Comprehensive Phase Diagrams of MoS Edge Sites Using Dispersion-Corrected DFT Free Energy Calculations Andrew S. Rosen, Justin M Notestein, and Randall Q. Snurr J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b02524 • Publication Date (Web): 08 Jun 2018 Downloaded from http://pubs.acs.org on June 12, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Comprehensive Phase Diagrams of MoS2 Edge Sites using Dispersion-Corrected DFT Free Energy Calculations Andrew S. Rosen, Justin M. Notestein*, Randall Q. Snurr* Department of Chemical and Biological Engineering, Northwestern University, Evanston, Illinois 60208, United States

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 38

ABSTRACT

A comprehensive set of surface phase diagrams addressing the catalytically relevant edges of the (100) surface of MoS2 catalysts is developed using dispersion-corrected density functional theory and ab initio thermodynamic modeling. The results of the temperature-dependent, free energybased thermodynamic model are presented over the full range of catalytically relevant temperatures and pressures, in addition to S- and H-coverages ranging from 0 – 100%. The results of this work allow for a full thermodynamic analysis to be performed at the conditions relevant to any promising reaction involving MoS2, ranging from hydrodesulfurization to dehydrogenation to electrocatalysis. Several methodological recommendations are discussed and implemented with the goal of improving the accuracy of the surface phase diagrams at minimal computational expense. A library of the most stable S and H adsorption modes is also developed so that linear scaling relationships can be used to correlate thermodynamic stability with kinetic activity. Applying the results to C-H bond activation of methane with an S2 oxidant, we predict S-coverages near 100% on the Mo- and S-edges to be thermodynamically favored and S monomers on edge sites with high S-coverages to be kinetically favorable. For H-abstraction on surface S atoms, the Mo-edge is also predicted to be more active than the S-edge.


2

Page 3 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION MoS2 materials have been employed for decades as industrial lubricants and as catalysts for the hydrodesulfurization (HDS) of crude oil and for sour water-gas shift reactions in the petrochemical industry,1,2 but in the past several years MoS2-based catalysts have demonstrated exceptional promise for a much larger number of applications. For instance, MoS2 can act as an effective catalyst for the hydrogen evolution reaction (HER),3,4 non-oxidative conversion of methane to ethylene,5 selective dehydrogenation of isobutane to isobutene,6 and dehydrogenative coupling of neat alcohols.7 For all of these X-H bond breaking/making reactions, it is essential to understand the structure and composition of the catalytic active sites under experimentallyrelevant operating conditions. This is critical for both the construction of an appropriate model system in computational studies and for a greater atomistic understanding of experimental data to enable data-driven catalyst design. Using density functional theory (DFT) and equilibrium thermodynamics, it is possible to develop an ab initio phase diagram of the MoS2 catalyst to identify the most favorable edge structures at a given set of reaction conditions. The first such thermodynamic model was developed by Raybaud et al.8 and only considered the reversible reaction of an S adatom with gas-phase H2. This model was improved by Bollinger et al.9 and Lauritsen et al.10 to account for the effect of adsorbed H atoms and later by Prodhomme et al.,11 who showed that the vibrational contributions of the edge adatoms can greatly impact the surface phase diagram. These models were developed in the context of specific operating conditions and pre-specified edge configurations relevant to HDS. The experimental conditions for other reactions involving an MoS2 catalyst can be quite diverse, spanning the large range of temperatures and pressures highlighted in Table 1. For reactions that are run at conditions different from these prior models,


3


Page 4 of 38

there is no clear way of predicting the most stable S- and H-coverages on the MoS2 edge sites based on the existing literature without ignoring potentially relevant surface coverages and excluding the vibrational contributions to the free energy. This naturally limits the ability to readily perform detailed, combined kinetic and thermodynamic investigations of novel reactions on MoS2 edge sites. In addition, much of the existing literature uses an “infinite stripe” model for the MoS2 catalyst. For reactions where multilayer MoS2 is expected to be catalytically relevant, the transferability of phase diagrams previously developed for monolayer MoS2 is not wellestablished. Table 1. Selected Reaction Conditions Involving MoS2 Catalysts Reaction

Temperature (K)

Hydrogen evolution reaction12

300

Dehydrogenative coupling of ethanol7

450 – 500

Hydrodesulfurization1

600 – 700

Alkane dehydrogenation6

600 – 900

Methane conversion using S2 as the oxidant5

1000 – 1300

Partial pressures (bar)

= 0.1 – 1

in equilibrium amountsa = 1 – 10

in equilibrium amountsa = 1 – 100

= 0.1 – 10

= 0.01 – 0.1b

in equilibrium amountsa / = 103 – 104

For reactions where H2S is not a reactant or product of the desired reaction, is expected to be present in small amounts from the reaction of gas-phase H2 with surface S atoms. This equilibrium value of is assumed to be small but is not readily measured. A prior HER study by Tsai et al., for instance, assumed a low range of 10-5 – 10-8 bar.12 a

b

Estimate based on the equilibrium conversion of 1 bar of isobutane using tabulated experimental data13 and assuming no side-reactions.


4

Page 5 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In this work, we develop a detailed and general ab initio thermodynamic model of the MoS2 edge sites that explicitly includes S- and H-coverages from 0% to 100% and that can be used for any temperature or pressure where the MoS2 catalyst is stable. In developing the thermodynamic model, we first identify the most stable configurations of S and H adatoms on the Mo- and Sedges at various coverages. With the most stable adsorption modes identified, we then determine the thermodynamically favorable edge structures at a wide range of temperatures and pressures. The data used to generate the phase diagrams, including DFT-computed electronic energies and all relevant vibrational frequencies, are reported so that the model can be easily extended to any set of experimental conditions not explicitly considered in this work. In this analysis, we also employ dispersion-corrected DFT, the first time such corrections have been considered for the MoS2 surface phase diagram. The inclusion of dispersion corrections affects the computed energetics and resulting thermodynamics, as MoS2 consists of alternatingly stacked monolayers that are stabilized by weak van der Waals (vdW) interactions that are not included in the standard Kohn-Sham formalism of DFT.14 Prior benchmark studies suggest that dispersion corrections greatly improve the structural properties of MoS2 and can be crucial in accurately predicting adsorption energies on periodic surfaces.15–18

METHODOLOGY Density Functional Theory Plane-wave DFT calculations were performed using the Vienna ab initio simulation package (VASP)19,20 v.5.4.1 with the electron exchange-correlation described by the Perdew-BurkeErnzerhof (PBE) functional.21 To account for vdW dispersion interactions, Grimme’s D3 dispersion correction with Becke-Johnson (BJ) damping was used for all calculations.22,23 The


5


Page 6 of 38

kinetic energy cutoff for the plane-wave basis set was chosen to be 650 eV. A Monkhorst-Pack24

4 × 2 × 2 -point grid was used to sample the Brillouin zone for all MoS2 edge structures, and

an 18 × 18 × 4 Monkhorst-Pack -point grid was used during the bulk MoS2 lattice parameter

optimization. Justifications for the energy cutoff and -point grid for modeling MoS2 are discussed in Tables S1 and S2 of the Supporting Information. Gaussian smearing of the band

occupancies with a smearing width of = 0.05 eV was implemented. The value of = 0.05 eV

was necessary to ensure that the total energies were converged in the extrapolated 0 K limit under the finite-temperature local-density approximation. The VASP 5.4 projector augmented-

wave (PAW) pseudopotentials25 were used for all calculations, where both the and semi-core

states were treated as valence for Mo, while the standard PAW pseudopotentials were used for all other elements. The accurate-precision keyword was specified in VASP to prevent

wraparound errors, and a dipole correction26 along the -axis was used to prevent spurious interactions between periodic images of the unit cell. Calculation details regarding other isolated and bulk molecules considered in this work can be found in the Supporting Information. Ionic relaxations were performed using a conjugate gradient algorithm until the net force on all individual atoms was less than 0.02 eV Å-1. A convergence criterion of 10-4 eV was set for the electronic energy minimization routine, except during the calculation of vibrational modes, in which a stricter convergence of 10-8 eV was set. Vibrational analyses were performed using central finite-difference approximations of the Hessian matrix with 0.01 Å displacements in each degree of freedom (DOF). The thermochemical calculations for all adsorbates were treated using the harmonic approximation,27 namely = +

'()

+

!" #$ %

−1

1


6

Page 7 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


'()

* = + ,

+ -

!" #$ %

− 1.

− ln -1 −

3 = − *

1

!" #$ % .2

3

2

where is the total DFT electronic energy, is the zero-point vibrational energy, and is

the vibrational energy of harmonic degree of freedom 4. Ideal gas statistical mechanics was used

for all gas-phase thermochemistry.28

The Atomic Simulation Environment (ASE)29 v.3.15.0 Python module was used to set up, run, and interpret all VASP calculations and to calculate the thermochemistry. Virtual NanoLab30 was used to generate the molecular visualizations, and the Matplotlib Python library31 was used to make the data visualizations. MoS2 Structural Parameters The most stable polytype of MoS2 exists in the space group P63/mmc (no. 194) and has three structural parameters to optimize within this constraint.32 Two of these are lattice constants: the Mo-Mo distance (5) and the interlayer separation (6/2). The third parameter is an internal degree

of freedom describing the distance between Mo and S in an MoS2 layer (7, given in units of 6).33

All three parameters were optimized simultaneously, and the details of the optimization

procedure are shown in Figures S1 and S2 of the Supporting Information. The minimum energy structural parameters were determined to be 5 = 3.15 Å, 6/2 = 6.07 Å, and 7 = 0.129. These

parameters are in excellent agreement with the experimental values of 5 = 3.16 Å, 6/2 = 6.15 Å,

and 7 = 0.129 from X-ray diffraction data.32 The interlayer separation only has a 1.5% deviation

from the experimental value, which is substantially better than the 10% error obtained using PBE without the D3(BJ) dispersion correction.15 MoS2 Model System


7


Page 8 of 38

The (100) surface of the hexagonal polytype of MoS2 is well-accepted as the catalytically active surface, with the basal plane assumed to be chemically inert for many reactions.3,4,12,34

Two chemically distinct edge types exist on the (100) surface: the (10180) edge and the (18010)

edge, referred to as the Mo-edge and S-edge, respectively. The S content on both edges is known to strongly depend on the reaction conditions, often varying from the crystallographic configuration.8,9,11,35–37 The crystallographic, as-cleaved (100) surface of bulk MoS2 used as the starting structure in

this work is shown in Figure 1 and consists of 72 atoms (24 Mo and 48 S atoms). The MoS2

layers are alternatingly stacked in the 9 direction, and the edge sites are oriented in the

direction for visual clarity. The simulation unit cell consists of 3 layers in the : direction, 1

bilayer in the 9 direction, and 4 layers in the direction, as adopted in prior computational

studies.8,38 During structural relaxations, the bottom two layers in the direction are fixed to mimic bulk constraints. To prevent artificial interactions with the neighboring periodic images,

10 Å of vacuum space is added above and below the MoS2 layers in the direction for a total of 20 Å of vacuum between each slab. The resulting unit cell has dimensions of 9.45 Å × 12.13 Å × 29.09 Å.


8

Page 9 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. The crystallographic (100) surface of MoS2 used as the simulation unit cell in this work. The Mo atoms are represented by the larger teal spheres, and the S atoms are represented by the smaller yellow spheres.

Ab initio Thermodynamic Model Grand Potential Formalism We consider a grand canonical ensemble wherein the number of particles can fluctuate and the

thermodynamically favorable edge structure is the one that minimizes the grand potential Φ at a

particular set of reaction conditions.39 In the case of an arbitrary MoSxHy edge structure, the grand potential can be evaluated as

Φ@AB C , E@A∗ , E∗ , E∗ = 3@AB C − G@A∗ E@A∗ − G∗ E∗ − G∗ E∗

4

where 3@AB C is the Helmholtz free energy of the edge structure under consideration, E ∗ is the chemical potential of an adsorbate 4, and G is the number of atoms of species 4. The value of

G@A∗ E@A∗ does not need to be explicitly computed for each configuration, as the number of Mo

atoms does not change between the structures considered in this work and so all relative grand

potentials will be independent of this term. For clarity, we also split the Helmholtz free energy

into two components: the DFT electronic energy @AB C and the vibrational contributions to the

HIJ Helmholtz free energy 3@A (including the zero-point vibrational energy). It can then be B C

stated that the grand potential of a given edge configuration is given by

HIJ − G∗ E∗ − G∗ E∗ Φ@AB C , E∗ , E∗ = @AB C + 3@A B C

5

In this analysis, the variations in surface coverages on the Mo- and S-edges are treated

independently, and the opposite edge is kept fixed at the crystallographic conditions: K = 0,


9


Page 10 of 38

K = 0 for the Mo-edge and K = 1, K = 0 for the S-edge. It is also assumed that each edge

structure is in the electronic ground state and that the translational and rotational degrees of freedom are negligible compared to the vibrational contributions. The configurational entropy is not expected to significantly impact the relative grand potentials and is therefore not considered. Vibrational Mode Analysis The vibrational modes of the adsorbates must be included in order to capture the significant thermal and entropic effects on the edge stability.11 This is in contrast with much of the prior literature that assumed 3@AB C ≈ @AB C in the expression for the grand potential or edge free

energy.9,10,36,37 Due to the lack of analytical frequency calculations in plane-wave DFT codes,

vibrational mode analyses are computationally expensive. Performing a detailed normal mode analysis including multiple surface atoms on each of the ~60 structures in this work is therefore undesirable. Hinneman et al.4 noted that the vibrational modes of H* on surface S atoms are independent of

both K and K . Under this approximation, the vibrational contribution to 3@AB C at a given temperature is a constant that can be scaled by the number of H atoms. We confirmed this finding and also found that this coverage independence holds for H* on surface Mo atoms as well, although there are slight differences depending on the local coordination environment. Based on a series of test calculations using the aforementioned DFT methodology, we define four sets of vibrational modes for H adatoms depending on the local environment: 1) H* on a S monomer directed toward the basal plane (595 cm-1, 645 cm-1, 2411 cm-1); 2) H* on a S monomer directed along the edge (504 cm-1, 603 cm-1, 2538 cm-1); 3) H* on a Mo-S-Mo interstitial site (757 cm-1, 1066 cm-1, 1243 cm-1); 4) H* on a Mo-Mo bridge site (550 cm-1, 798


10

Page 11 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cm-1, 1116 cm-1). These four modes are visually depicted in Figure 2. We further take advantage of this simplification by hypothesizing that the change in K negligibly impacts the S*

vibrational modes so that the vibrational frequencies of S* only need to be calculated for K = 0.

The vibrational frequencies of S* for K = 0 − 1 and K = 0 can be found in the Supporting

Information. The effectively decoupled vibrational modes are likely due to the large mass differences between H, S, and Mo.

Figure 2. Schematic showing four common adsorption modes of H* on surface Mo and S atoms. To summarize, with this set of vibrational frequencies, the vibrational contribution to the

HIJ Helmholtz free energy for a given edge structure, 3@A , can be calculated as B C HIJ HIJ = 3HIJ 3@A ∗ + G∗ 3∗ B C

6

where 3HIJ ∗ is the vibrational contribution to the Helmholtz free energy of the S adatoms at a

given K (with K = 0), and 3HIJ ∗ is the coverage-independent vibrational contribution to the

Helmholtz free energy of a single H adatom. The accuracy of this decoupled, scaled approach for

calculating the vibrational contributions to the free energy compared to a full vibrational analysis is reported in Table S3 of the Supporting Information for representative structures. The error is below that expected from DFT and the harmonic approximation. Expressions for NO∗ and NP∗


11


Page 12 of 38

When the MoS2 catalyst is in the presence of H2, the S atoms on the catalytic edge sites can reversibly react to form H2S via

S ∗ + HS g ↔ HS S g + ∗

7

where * denotes a S-vacancy. At thermodynamic equilibrium, the chemical potentials can be equated such that

E ∗ = E − E

8

1 H g +∗ 2 S

9

Similarly, adsorbed H atoms can reversibly form gas-phase H2 via H∗ ↔

where * denotes an H-adsorption site. Under equilibrium conditions, E ∗ =

1 E 2

10

It is possible to write an analogous expression for E∗ and E∗ in terms of other species that are

more suitable for a given reaction. In the case of isobutane dehydrogenation to isobutene,6 for instance, it may be more natural to think of E∗ as a function of EWX YZ and EWX [ , whereas for

methane coupling using gaseous S2 as the oxidant,5 E∗ could instead be written as a function of E . For consistency, we present the results in terms of E and E . The chemical potential of a gas-phase species 4 is given by ] E = E∘ + + ln - ∘ . ]

11

where E∘ is the reference chemical potential of the gas, ] is the fugacity of the gaseous species,

and ] ∘ is the reference fugacity (chosen to be 1 bar). At low to moderate pressures where the

gas-phase species can be assumed to be ideal, the gas-phase fugacity approaches the partial

pressure . We present all phase diagrams in terms of fugacity for generality, but the fugacity

can be assumed to be approximately equal to for all the reaction conditions listed in Table 1. ACS Paragon Plus Environment

12

Page 13 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Even for high-pressure HDS conditions, the fugacity coefficient for 100 bar of H2 at 650 K is 1.02 as calculated from the Peng-Robinson equation of state, indicating nearly ideal behavior.

With these thermodynamic relations, E∗ and E∗ can be written as the following under

equilibrium conditions:

∘ ∘ E∗ = E − E + + ln ^

E ∗ =

] _ ]

] 1 ∘ `E + + ln ^ ∘ _a 2 ]

The reference chemical potentials E∘ were calculated from9

12

13

E∘ = b∘ − b∘ 0 K + + − *∘

14

where b∘ is the standard molar enthalpy, is the zero-point energy, is the energy obtained

from DFT, and *∘ is the standard molar entropy. The extended Shomate equation with

experimental parameters from the NIST chemistry webbook40 were used to express b∘ and *∘ as a function of , and the standard state enthalpies at 0 K were obtained from the NIST-JANAF

thermochemical tables.41 Vibrational frequencies used to calculate were obtained from the

CRC Handbook.42

In prior literature on the MoS2 surface phase diagram, the absolute chemical potential E∗ is

written with the DFT electronic energy (per S atom) of the orthorhombic d-phase of bulk S used

as a reference. While reference states are commonly the most stable pure allotrope at 1 bar and

298 K, this choice is arbitrary. For the sake of DFT calculations, Jefg can be a problematic

reference state for E∗ , as d-S consists of S8 rings stabilized by vdW interactions, which are

known to not be properly described by DFT. Any error introduced in the choice of reference state will exponentially impact the accuracy of calculating the fugacities from equations (12) and


13


Page 14 of 38

(13) that correspond to a given chemical potential. We instead use H2S and H2 as the reference state such that

15

ΔE∗ = E∗ − i − j

which will generally be more accurate to calculate regardless of the choice of density functional. The reference state for calculating ΔE∗ is chosen to be gas-phase H2 as in prior studies such that 1 ΔE∗ = E∗ − 2

16

MoS2 is only stable within a certain range of chemical potentials. The lower and upper limits

of ΔE∗ represent the transition to forming bulk Mo and bulk S, respectively. Using the same

methodology described previously,9 this corresponds to approximately −1.26 eV ≤ ΔE∗ ≤

0.13 eV. Similarly, the upper-bound on ΔE∗ is approximately 0 eV, corresponding to the point where the formation of gas-phase H2 is preferable over H*. The lower-bound for ΔE∗ is

configuration-dependent and corresponds to the limit where H atoms will adsorb on the MoS2

edge.9 For data presentation purposes, we set the lower-bound to be the ΔE∗ value that corresponds to ] /] ∘ = 10-16, representative of ultrahigh vacuum (UHV) conditions.

RESULTS AND DISCUSSION Stable Adsorption Modes

We first determine the preferential S* adsorption sites at various K values while no H atoms

are on the surface. Using the most stable S* configurations at each K value, we then consider the addition of H atoms and determine their preferential adsorption sites. This stepwise procedure is

done to ensure that the structures are at – or at least close to – the global minimum of the potential energy surface.


14

Page 15 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The most stable S* configurations at different values of K were computationally investigated

by Raybaud et al.8 nearly two decades ago using the PW91 functional43 and ultrasoft pseudopotentials. To test whether the same results are found using the computational details in this work, the relative stabilities for different S* configurations on both the Mo- and S-edges were reinvestigated using the set of proposed S-adsorption modes by Raybaud et al.8 as well as

those investigated by Prodhomme et al.11 for 100% S-coverage on the S-edge. The value of K on

a given edge is defined as

K =

G 2G@A

17

where G represents the number of terminal S atoms and G@A represents the number of Mo atoms

on the edge (of which there are 3 in the simulation unit cell).

With the most stable S* configurations known, we consider the sequential addition of H atoms

and determine the most stable adsorption modes at varying K and K . The value of K is defined as

K =

G G@A

18

where G is the number of adsorbed H atoms, as previously adopted by Tsai et al.12 We considered the set of possible adsorption sites identified by Kronberg et al.44 as well as select

additional configurations on the S-edge. The definition of K was chosen so that it is with respect to the same quantity for each edge structure (i.e. G@A = 3). As a result, some edge structures (e.g.

the S-edge with K = 1) have S adatoms that are not fully saturated by H atoms at K = 1.

However, these high H-coverages are not expected to be thermodynamically favorable within the stable range of ΔE∗ .


15


Page 16 of 38

The most stable adsorption sites and the corresponding DFT electronic energies are shown in Table 2 for the Mo-edge and Table 3 for the S-edge (the energy and coordinates of ~150 tested structures are included in the Supporting Information). The most stable S* configurations are in agreement with the previously calculated results of Raybaud et al.8 Focusing on the S-edge with

K = 1, Prodhomme et al.11 predict that an alternating sequence of S monomers/dimers is 0.99 eV

lower in energy than the full monomer configuration, at least for monolayer MoS2. As shown in Table S4 of the Supporting Information, we instead find the full monomer configuration to be more stable by 0.17 eV for bulk MoS2. It is also noted that due to the periodicity of the

simulation unit cell, the S-edge with K = 0.5 does not have a complete zig-zag configuration as noted by Schweiger et al.,45 although this is not expected to influence the results, as justified in Table S5 of the Supporting Information. Table 2. Stable S* and H* Adsorption Modes on the Mo-edge and Corresponding DFT Electronic Energies as a Function of K and K a

K = 0 K = 0.17 K = 0.33

K = 0

K = 0.33

K = 0.67

K = 1

Δ = -0.89

Δ = -0.78

Δ = -0.41

Δ = -2.91

Δ = -0.60

Δ = -0.20

Δ = 0.03

Δ = -2.37

Δ = -0.19

Δ = 0.15

Δ = 0.12


16

Page 17 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


K = 0.5 K = 0.67 K = 0.83 K = 1

Δ = -0.93

Δ = -0.64

Δ = 0.22

Δ = 0.17

Δ = -0.21

Δ = -0.77

Δ = 0.14

Δ = 0.49

Δ = -0.08

Δ = 0.03

Δ = 0.02

Δ = -0.17

Δ = -0.36

Δ = 0.34

Δ = 0.15

Δ = 0.43

Δ is the differential S-adsorption energy (in eV) to go from MoSx-1Hy to MoSxHy, and Δ is the differential H-adsorption energy (in eV) to go from MoSxHy-1 to MoSxHy a

Table 3. Stable S* and H* Adsorption Modes on the S-edge and Corresponding DFT Electronic Energies as a Function of K and K a

K = 0 K = 0.17 K = 0.33

K = 0

K = 0.33

K = 0.67

K = 1

Δ = -1.55

Δ = -1.63

Δ = -0.71

Δ = -3.84

Δ = -1.06

Δ = -1.00

Δ = 0.09

Δ = -3.61

Δ = -0.86

Δ = -0.17

Δ = 0.13


17


K = 0.5 K = 0.67 K = 0.83 K = 1

Page 18 of 38

Δ = -2.31

Δ = -0.28

Δ = -0.36

Δ = -0.34

Δ = -0.98

Δ = -0.55

Δ = -0.42

Δ = -0.50

Δ = -1.05

Δ = -0.41

Δ = -0.65

Δ = -0.59

Δ = -0.94

Δ = -0.73

Δ = -0.57

Δ = -0.51

Δ is the differential S-adsorption energy (in eV) to go from MoSx-1Hy to MoSxHy, and Δ is the differential H-adsorption energy (in eV) to go from MoSxHy-1 to MoSxHy a

In terms of H-adsorption, H atoms are predicted to adsorb to bare Mo sites before adsorbing to terminal S atoms on the surface of the Mo-edge. We also find that it is energetically unfavorable for H atoms to adsorb and break an S-S dimer present on the Mo-edge unless it is the only

remaining adsorption site. Similarly, for the case of K = 1 on the Mo-edge, the lowest energy H-adsorption mode is on the Mo-Mo bridge sites as opposed to the dimerized S-S bridge sites. It

has been previously reported that it is energetically favorable to break the S-S dimer for K = 1

on the Mo-edge only after significant reconstruction of the surface adatoms,44 although the barrier for this bond cleavage is greater than 0.8 eV.12 As a result, H-adsorption on S-S dimers followed by significant surface reconstruction was not considered in the present study for the Mo-edge with K = 1.


18

Page 19 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


On the S-edge, it is also found that H atoms tend to adsorb to bare Mo sites first. One

exception is for K = 0.67 at K = 0.33, although the energy difference between the structure

shown in Table 3 and the structure with H* on an interstitial Mo site is only 0.11 eV. On the Sedge with K > 0.5, H atoms on surface S atoms preferentially orient toward the adjacent basal

plane rather than an on-top adsorption mode, in agreement with prior studies.11,12 The favorable adsorption modes on the Mo-edge are generally consistent with those identified by Kronberg et al.,44 while most of the favorable structures on the S-edge shown in Table 3 were not investigated in this prior work.

The differential S-adsorption energies (at K = 0) defined as

Δ = @AB C + − @ABnY C −

19

and the differential H-adsorption energies defined as

1 Δ = @AB C − @AB CnY − 2

20

are also reported in Table 2 and Table 3. Using these definitions, Δ at K = 0.33, for instance,

refers to the energy to go from K = 0.17 to K = 0.33 (an analogous interpretation holds for the

definition of Δ ). For the H-free surface, we find that S-adsorption is always energetically

favorable on both the Mo- and S-edges. Unlike the findings of Raybaud et al.,8 S-adsorption

between K = 0.5 − 0.83 on the Mo-edge is not found to be endothermic, although the process is

less energetically favorable than at lower S-coverages. This trend of continuously favorable Sadsorption is attributed to the inclusion of dispersion corrections. The electronic energies of the geometry-optimized structures using the PBE functional without the D3(BJ) correction are analogous to those reported by Raybaud et al. (in this case, endothermic adsorption for K =

0.67 − 1). As shown in Figure S3 of the Supporting Information, using the dDsC dispersion-

correction scheme46 instead of the D3(BJ) scheme does not change the reported trends. In terms


19


Page 20 of 38

of H-adsorption, we find that as the H-coverage increases, H-adsorption tends to become weaker for many of the surface coverages on the Mo- and S-edges. No clear overarching trends are observed as a function of K .

With this extensive library of MoS2 edge structures and corresponding energies, the phase

diagram for the MoS2 edges can be generated by calculating Φ@AB C over a range of E∗ and

E ∗ .

MoS2 Edge Phase Diagrams We first consider the MoS2 edge phase diagram at 300 K, as shown in Figure 3. For brevity,

we use the abbreviation “K /K ” (e.g. Mo-50/33 for K = 0.5, K = 0.33 on the Mo-edge) to denote a given edge coverage. On the Mo-edge, we find that even at low ΔE∗ values, the surface does not exist in the crystallographic K = 0 state. Over the range of physically accessible ΔE∗

values, the Mo-edge can exist with K = 0.33 − 1 and tends to have low H-coverages of

K = 0 − 0.33. On the S-edge, much of the accessible portion of the phase diagram results in fully saturated surface Mo atoms with K = 1 and an H-coverage that is strongly dependent the

pressure of H2. Only at very low values of ΔE∗ , far below most kinetic studies involving an MoS2 catalyst, can the S-edge exist with K < 1. We also investigated K = 1 on the S-edge with higher H-coverages than those shown in Table 3 and found that K > 1 is not

thermodynamically favorable for the stable values of ΔE∗ .


20

Page 21 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Calculated MoS2 surface phase diagrams showing the most thermodynamically stable

edge configurations at 300 K as a function of the chemical potential of adsorbed sulfur ΔE∗ ,

chemical potential of adsorbed hydrogen ΔE∗ , and fugacities ] of gaseous H2S and H2. The reference fugacity ] ∘ is 1 bar. (a) Mo-edge. (b) S-edge.

HER studies and many forms of characterization, such as scanning tunneling microscopy (STM), X-ray photoelectron spectroscopy, and Fourier-transform infrared spectroscopy, are commonly done around 300 K, so it is worth briefly considering the most stable edge structures

for these conditions. For HER, the definitions of E∗ and E∗ can be readily modified to include an additive term for the applied potential based on the computational hydrogen electrode model if needed.47 For simplicity, we consider the case of low overpotential where this factor does not need to be included.12 In the operating regime for HER (refer to Table 1), the Mo-edge is most stable as Mo-50/33, while the S-edge is most stable as S-100/100. This is essentially the same as what Tsai et al. previously proposed.12 During STM imaging (and many other characterization techniques), the catalyst is studied under UHV. Based on the thermodynamic analysis, it is expected that the catalyst under UHV conditions would exist with full S-coverage on the S-edge


21


Page 22 of 38

and either full S-coverage or 50% S-coverage on the Mo-edge depending on the ratio of H2S to H2. This is in direct agreement with experimental STM images of hexagonal MoS2 monolayers.10 By considering a larger range of surface coverages than prior STM-focused studies for HDS catalysis,9,10 the present model suggests a few additional yet notable regions of stability at 300 K. When there is a low ratio of H2S to H2, sub-50% S-coverage is likely to be stable on the Moedge. At H2 partial pressures slightly higher than those used in STM imaging, we also predict

that a combination of S monomers and S-S dimers (i.e. K = 0.67) can be present on the Moedge. The stable S-coverages on the S-edge are similar to those proposed in prior models,9,10

although this is the first time that the theoretical phase transitions between low and high Hcoverages have been identified. The phase diagrams at 650 K for the Mo- and S-edges are shown in Figure 4. On the Mo-edge,

it is possible to have an edge configuration with K = 0.5 and K = 0 without resorting to UHV

conditions, in contrast with the 300 K phase diagram. At common pressures for HDS, we find

that the most likely Mo-edge configuration is Mo-50/33 or Mo-67/33 depending on the experimental conditions. While this goes against much of the older literature that suggested K = 1 is most stable,8,34 it is in reasonable agreement with the results of Prodhomme et al. who

suggest that Mo-50/50 is favorable.11 The model of Prodhomme et al. only accounted for

K = 0.37 − 0.5 on the Mo-edge,11 so the present model also uniquely highlights the large region of stability for the Mo-100/0 edge when the partial pressure of H2S is at least as high as

the partial pressure of H2. On the S-edge, we predict that the S-100/100 edge is most stable under HDS conditions. This is in agreement with the model developed by Bollinger et al.9 but less so with other studies that suggest K = 0.5 is likely.8,11 When calculating Jefg , we found that

the energetic contribution from the D3(BJ) dispersion correction was approximately -0.2


22

Page 23 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


eV/atom with the PBE functional, which changes the corresponding ] /] ratio by 2 orders of

magnitude at 650 K. Incorporating this correction in prior models (or, preferably, using −

as the energy reference for E∗ ) almost entirely accounts for the underestimated K values on

the S-edge in prior studies of HDS conditions.




The phase diagram at 1000 K is shown in Figure 5. At such a high temperature, the Mo-33/0

edge effectively becomes unstable, and the Mo-edge exists with either K = 0.5 or K = 1 for most of the phase diagram. Similarly, the intermediate S-coverages on the S-edge become less

likely under equilibrium conditions, with K = 1 taking up most of the accessible phase diagram.

Even without knowing the exact partial pressure of H2S and H2 in the system, the phase diagrams shown in Figure 5 significantly narrow down the possible S-coverages on both the Mo- and Sedges. To summarize the results from 300 K to 1000 K, the thermodynamically favorable surface coverages over the wide range of reactions conditions listed in Table 1 are presented in Table 4.


23


Page 24 of 38




Table 4. Predicted Stable K and K on the Mo-edge and S-edge for the Range of Experimental Reaction Conditions Outlined in Table 1. Reaction Hydrogen evolution reaction

Mo-edge (K /K ) 50/33

S-edge (K /K ) 100/100

Dehydrogenative coupling of 50/33 ethanol

100/100

Hydrodesulfurization

50/33, 67/33

100/100

Alkane dehydrogenation

33/0, 50/0, 50/33

50/0, 100/33, 100/67

Methane conversion using S2 100/0 as the oxidant

100/0

Considering the phase diagrams over this large range of temperatures, a number of generalizations can be made. Increasing the temperature has two major consequences. The first is


24

Page 25 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a direct result of the definitions of the chemical potentials: as increases, the range of fugacities

corresponding to a given range of chemical potentials decreases. At 300 K, a range of ΔE∗

values spanning 0.4 eV corresponds to a range of ] values spanning nearly 12 orders of

magnitude. At 1000 K, this same 0.4 eV range corresponds to just 4 orders of magnitude in ] .

Because of this, at high temperatures the MoS2 edge sites have 0% H-coverage over a much HIJ larger range of ] . The second consequence of changing is the effect of 3@A on Φ@AB q B C

and thereby the boundaries of the phase diagram at a given range of chemical potentials. Increasing the temperature from 300 K to 1000 K causes the phase boundaries to shift by as much as -0.4 eV in ΔE∗ . A similar phenomenon can be inferred from previously published

HDS-relevant free energy diagrams.11 This shift in the phase boundaries is what results in the inaccessibility of the Mo-33/0 and S-50/0 regions at high temperatures, greatly increasing the likelihood that both edges have high S-coverages. There is less of an impact on the phase

boundaries in the ΔE∗ dimension: going from 300 K to 1000 K generally results in a -0.1 eV shift in the phase boundaries, although the phase transition between Mo-67/33 and Mo-100/0

shifts by approximately 0.2 eV. Interestingly, changing the temperature does not seem to introduce new stable edge configurations. For instance, the Mo-83/0 and S-67/0 edges do not appear to be thermodynamically favorable at any temperature. In addition, the Mo-edge cannot

exist with K ≤ 0.17, and the S-edge cannot exist with K ≤ 0.33 regardless of temperature.

There are no upper-limits on K for either edge, and K can take on any value between 0% and

100% depending on the operating conditions.

The applicability of the thermodynamic data extends beyond predicting the most stable surface structures for the aforementioned reactions. Any set of temperatures and pressures can easily be investigated using this thermodynamic model, as the Supporting Information contains all


25


Page 26 of 38

relevant energies and vibrational frequencies along with a Python script to quickly perform the thermodynamic calculations. For convenience, phase diagrams for the Mo- and S-edges from 300 K – 1000 K in increments of 100 K are also reported in Figure S4 of the Supporting Information. Structure-Property Relationships In addition to providing insight into the thermodynamic stability of different edge structures, the free energies reported in this work can be combined with known catalytic structure-property relationships to provide significant insight into the reaction kinetics. For instance, it has been previously reported that Brønsted-Evans-Polanyi relationships exist for H2 and H2S associative desorption,11 that HER activity is correlated with the Gibbs free energy of H-adsorption via the Sabatier principle,48 and that the barrier for activating the C-H bond of methane is inversely related to the binding strength of surface S atoms on metal sulfides.5 From these examples, we consider the catalytic conversion of methane to ethylene in the presence of gaseous S2 to highlight how the data in the present study can be used to quickly predict which edge structures may be catalytically promising for the conversion of methane (and potentially other reactions involving a similar mechanism for H-abstraction). In the combined experimental and theoretical work of Zhu et al.,5 it was shown that MoS2 can activate the C-H bond of methane, forming a surface-stabilized CH3 species. The CH3 intermediate can go through further H-abstractions, and it was hypothesized that the coupling of nearby CH2 radicals on MoS2 edge sites can result in the production of ethylene. Among the four metal sulfides tested for this reaction, MoS2 had the highest methane conversion but the lowest selectivity to ethylene. The authors found that the M-S bond strength (M = Mo, Pd, Ru, Ti) is directly correlated to the barrier for C-H activation of methane.


26

Page 27 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The M-S bond strength is equivalent to the differential S-desorption energy and is calculated as the energy lost when a surface S atom reacts with H2 and leaves as H2S. Focusing on MoS2, the

standard Gibbs free energy at 1000 K for this process is shown in Figure 6 as a function of K .

The Mo-S bond energies are based on the most stable H-free edge sites shown in Table 2 and Table 3. Additional Mo-S bond energies are reported in Table S6 of the Supporting Information.

Overall, similar trends are observed on both the Mo- and S-edges: as S atoms are lost as H2S from the MoS2 surface, the Mo-S bond energy tends to increase. Based on the Mo-S bond energy descriptor and the energies reported in Figure 6, it is likely that the S adatoms on edge sites with higher S-coverages are among the most active for breaking the strong C-H bond of methane. From the equilibrium phase diagram at 1000 K shown in Figure 5, K = 1 is the most thermodynamically stable S-coverage on the Mo- and S-edges at the high ΔE∗ values associated

with using S2 as an oxidant. This implies that the most thermodynamically stable edge structures may also have S adatoms that are the most active for H-abstraction, an inherent requirement for the observed conversion of methane. One possible exception is the case of S-S dimers on the surface of the Mo-edge at high S-coverages, as an H-abstraction mechanism requiring the simultaneous cleavage of an S-S bond is unlikely to result in a low activation energy.


27


Page 28 of 38

Figure 6. Mo-S bond energy as a function of K for the S-edge (red lines) and Mo-edge (blue lines). The value at a given K refers to the energy difference (dashed line: electronic energy,

solid line: standard Gibbs free energy at 1000 K) associated with removing a surface S atom as H2S at that S-coverage. In the context of oxidative coupling of methane to ethylene, too low an Mo-S bond strength results in large barriers for CH2 coupling as a consequence of the Sabatier principle.5 Decreasing the S-coverage on either edge to below half a monolayer by decreasing the partial pressure of S2 is hypothesized here to increase the otherwise poor selectivity to ethylene that was experimentally observed with an MoS2 catalyst, albeit with a significant decrease in methane conversion.5 By the same logic, we predict that S adatoms on the Mo-edge are generally more active for the C-H bond activation step, whereas the S-edge is more active for CH2 coupling to ethylene. It then becomes clear that there are significant structure-property relationships that can be used to finely tune the conversion and selectivity by modifying the local surface coverage and edge type. This type of combined thermodynamic and kinetic analysis can be readily performed for other reactions as well, so long as linear scaling relationships are known.48 CONCLUSIONS We have developed a comprehensive set of surface phase diagrams of the Mo- and S-edges on the (100) surface of MoS2 using density functional theory with dispersion corrections. We first determined the most stable S* and H* adsorption modes and their corresponding electronic

energies for K and K from 0 – 100%. The vibrational modes of S* and H* were then

considered in the development of the thermodynamic model to improve the accuracy of the predicted phase boundaries as a function of temperature.


28

Page 29 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Using the stable structures, electronic energies, and adsorbate vibrational modes, surface phase diagrams were developed over the full range of relevant operating conditions for several important catalytic reactions. We found that increasing the temperature does not result in the emergence of new stable coverages, although it does shift the phase boundaries by as much as 0.4 eV in ΔE∗ . This reduces the thermodynamic stability of lower S-coverages at higher

temperatures, making the K ≤ 0.33 region on the Mo-edge and K ≤ 0.83 region on the S-edge

thermodynamically unfavorable at 1000 K. In terms of adsorbed H atoms, the Mo-edge is most

stable with either K = 0 or K = 0.33 depending on ΔE∗ , whereas the S-edge can exhibit a

wider range of H-coverages depending on ΔE∗ . This is true across the full range of temperatures

considered in this work, as the phase boundaries generally shift by -0.1 eV in the ΔE∗

dimension when the temperature is increased from 300 K to 1000 K.

In constructing the ab initio thermodynamic phase diagrams, the inclusion of vibrational degrees of freedom was found to be necessary for capturing the boundaries for each of the phase

transitions, particularly as a function of ΔE∗ . In materials like MoS2 where the large mass differences between H, S, and Mo result in nearly decoupled vibrational modes of H adatoms

and their S or Mo adsorption sites, only a few additional calculations are required. This is especially true if the vibrational frequencies associated with H-adsorption are independent of surface coverage, as we confirmed for MoS2. We considered a few operating conditions as examples to highlight the utility of the phase diagrams presented in this work. For HER at 300 K with low overpotential, the Mo-50/33 and S100/100 edges are most stable. For UHV conditions at 300 K (e.g. STM imaging), the Scoverage is either 50% (S monomers) or 100% (S-S dimers) on the Mo-edge and 100% on the Sedge. For HDS, the most stable coverages were found to be Mo-50/33 or Mo-67/33 (depending


29


Page 30 of 38

on the pressure of H2) and S-100/100. Finally, applying the results to an example involving MoS2 for the conversion of methane to ethylene, we predict that 100% S-coverage is the most stable for both the Mo- and S-edges and that S-monomers at high S-coverages are optimal for the conversion of methane when using Mo-S bond strength as a descriptor for the C-H activation barrier. In addition, for H-abstraction via surface S atoms, the Mo-edge is expected to be more active than the S-edge. With the data presented in this study, the MoS2 surface stability analysis is extended beyond the commonly studied conditions and edge configurations of the HDS catalyst so that it can be more generally applied to any reaction involving surface S* and H* atoms. The comprehensive phase diagrams and tabulated energies presented in this work open up the possibility to readily perform a full thermodynamic analysis for new reactions involving an MoS2 catalyst and can provide preliminary insight into structure-property relationships relevant to experimental trends in catalytic activity. ASSOCIATED CONTENT Supporting Information Additional details for DFT calculations, justification of vibrational mode assumptions, discussion of unit cell periodicity, additional Mo-S bond energies, summarized methodology, phase diagrams from 300 K – 1000 K in increments of 100 K (PDF); Python script to reproduce the stability analysis in this work at user-specified conditions, an Excel sheet containing the DFT electronic energies and vibrational modes, atomic coordinates of the optimized structures considered in this work (.zip). AUTHOR INFORMATION


30

Page 31 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Corresponding Author E-mail: [email protected]. Tel.: 847-467-2977 E-mail: [email protected]. Tel.: 847-491-5357 Author Contributions The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. ACKNOWLEDGMENTS This work was supported as part of the Inorganometallic Catalyst Design Center, an Energy Frontier Research Center funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), under Award DE-SC0012702. A.S.R. acknowledges Government support under contract FA9550-11-C-0028 and awarded by the Department of Defense, Air Force Office of Scientific Research, National Defense Science and Engineering Graduate (NDSEG) Fellowship, 32 CFR 168a. The authors acknowledge computing support through the resources and staff contributions provided for the Quest high-performance computing facility at Northwestern University as well as the resources provided by the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. REFERENCES (1)

Topsøe, H.; Clausen, B. S.; Massoth, F. E. Hydrotreating Catalysis. In Catalysis; Springer, 1996; pp 1–269.

(2)

Łaniecki, M.; Małecka-Grycz, M.; Domka, F. Water-Gas Shift Reaction over Sulfided


31


Page 32 of 38

Molybdenum Catalysts: I. Alumina, Titania and Zirconia-Supported Catalysts. Appl. Catal. A Gen. 2000, 196, 293–303. (3)

Jaramillo, T. F.; Jørgensen, K. P.; Bonde, J.; Nielsen, J. H.; Horch, S.; Chorkendorff, I. Identification of Active Edge Sites for Electrochemical H2 Evolution from MoS2 Nanocatalysts. Science. 2007, 317, 100–102.

(4)

Hinnemann, B.; Moses, P. G.; Bonde, J.; Jørgensen, K. P.; Nielsen, J. H.; Horch, S.; Chorkendorff, I.; Nørskov, J. K. Biomimetic Hydrogen Evolution: MoS2 Nanoparticles as Catalyst for Hydrogen Evolution. J. Am. Chem. Soc. 2005, 127, 5308–5309.

(5)

Zhu, Q.; Wegener, S. L.; Xie, C.; Uche, O.; Neurock, M.; Marks, T. J. Sulfur as a Selective “Soft” Oxidant for Catalytic Methane Conversion Probed by Experiment and Theory. Nat. Chem. 2013, 5, 104–109.

(6)

Wang, G.; Li, C.; Shan, H. Highly Efficient Metal Sulfide Catalysts for Selective Dehydrogenation of Isobutane to Isobutene. ACS Catal. 2014, 4, 1139–1143.

(7)

McCullough, L. R.; Childers, D. J.; Watson, R. A.; Kilos, B. A.; Barton, D. G.; Weitz, E.; Kung, H. H.; Notestein, J. M. Acceptorless Dehydrogenative Coupling of Neat Alcohols Using Group VI Sulfide Catalysts. ACS Sustain. Chem. Eng. 2017, 5, 4890–4896.

(8)

Raybaud, P.; Hafner, J.; Kresse, G.; Kasztelan, S.; Toulhoat, H. Ab Initio Study of the H2– H2S/MoS2 Gas–Solid Interface: The Nature of the Catalytically Active Sites. J. Catal. 2000, 189, 129–146.

(9)

Bollinger, M. V.; Jacobsen, K. W.; Nørskov, J. K. Atomic and Electronic Structure of MoS2 Nanoparticles. Phys. Rev. B 2003, 67, 85410.


32

Page 33 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(10)

Lauritsen, J. V.; Bollinger, M. V.; Lægsgaard, E.; Jacobsen, K. W.; Nørskov, J. K.; Clausen, B. S.; Topsøe, H.; Besenbacher, F. Atomic-Scale Insight into Structure and Morphology Changes of MoS2 Nanoclusters in Hydrotreating Catalysts. J. Catal. 2004, 221, 510–522.

(11)

Prodhomme, P. Y.; Raybaud, P.; Toulhoat, H. Free-Energy Profiles along Reduction Pathways of MoS2 M-Edge and S-Edge by Dihydrogen: A First-Principles Study. J. Catal. 2011, 280, 178–195.

(12)

Tsai, C.; Chan, K.; Abild-Pedersen, F.; Nørskov, J. K. Active Edge Sites in MoSe2 and WSe2 Catalysts for the Hydrogen Evolution Reaction: A Density Functional Study. Phys. Chem. Chem. Phys. 2014, 16, 13156–13164.

(13)

Dean, J. A. Lange’s Handbook of Chemistry, 15th ed.; McGraw-Hill Professional: New York, 1998.

(14)

Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133.

(15)

Peelaers, H.; Van de Walle, C. G. First-Principles Study of van Der Waals Interactions in MoS2 and MoO3. J. Phys. Condens. Matter 2014, 26, 305502.

(16)

Björkman, T.; Gulans, A.; Krasheninnikov, A. V; Nieminen, R. M. Are We van Der Waals Ready? J. Phys. Condens. Matter 2012, 24, 424218.

(17)

Kerber, T.; Sierka, M.; Sauer, J. Application of Semiempirical Long-Range Dispersion Corrections to Periodic Systems in Density Functional Theory. J. Comput. Chem. 2008, 29, 2088–2097.


33


(18)

Page 34 of 38

Ganji, M. D.; Hosseini-Khah, S. M.; Amini-Tabar, Z. Theoretical Insight into Hydrogen Adsorption onto Graphene: A First-Principles B3LYP-D3 Study. Phys. Chem. Chem. Phys. 2015, 17, 2504–2511.

(19)

Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169–11186.

(20)

Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector AugmentedWave Method. Phys. Rev. B 1999, 59, 1758–1775.

(21)

Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868.

(22)

Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104.

(23)

Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456–1465.

(24)

Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13, 5188–5192.

(25)

Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953–17979.

(26)

Makov, G.; Payne, M. C. Periodic Boundary Conditions in Ab Initio Calculations. Phys. Rev. B 1995, 51, 4014–4022.

(27)

McQuarrie, D. A. Statistical Mechanics; University Science Books: Sausalito, 2000.


34

Page 35 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(28)

Cramer, C. Essentials of Computational Chemistry, 2nd ed.; Johny Wiley and Sons, Ltd.: West Sussex, 2004.

(29)

Larsen, A.; Mortensen, J.; Blomqvist, J.; Castelli, I.; Christensen, R.; Dulak, M.; Friis, J.; Groves, M.; Hammer, B.; Hargus, C.; et al. The Atomic Simulation Environment—A Python Library for Working with Atoms. J. Phys. Condens. Matter 2017, 29, 273002.

(30)

Virtual NanoLab version 2017.1, QuantumWise A/S www.quantumwise.com.

(31)

Hunter, J. D. Matplotlib: A 2D Graphics Environment. Comput. Sci. Eng. 2007, 9, 90–95.

(32)

Bronsema, K. D.; De Boer, J. L.; Jellinek, F. On the Structure of Molybdenum Diselenide and Disulfide. ZAAC

(33)

J. Inorg. Gen. Chem. 1986, 540, 15–17.

Grau-Crespo, R.; López-Cordero, R. Comment on “Ab Initio Study of MoS2 and Li Adsorbed on the (1010) Face of MoS2” by V. Alexiev, R. Prins and Th. Weber, Phys. Chem. Chem. Phys., 2000, 2, 1815, and “DFT Study of MoS2 and Hydrogen Adsorbed on the (1010) Face of MoS2” by V. Alexiev, R. Prins. Phys. Chem. Chem. Phys. 2002, 4, 4078–4079.

(34)

Byskov, L. S.; Nørskov, J. K.; Clausen, B. S.; Topsøe, H. DFT Calculations of Unpromoted and Promoted MoS2-Based Hydrodesulfurization Catalysts. J. Catal. 1999, 187, 109–122.

(35)

Plazenet, G.; Cristol, S.; Paul, J.-F.; Payen, E.; Lynch, J. Sulfur Coverage and Structural Disorder in γ-Alumina-Supported MoS2 Crystallites: An in Situ EXAFS Study. Phys. Chem. Chem. Phys. 2001, 3, 246–251.

(36)

Cristol, S.; Paul, J. F.; Payen, E.; Bougeard, D.; Clémendot, S.; Hutschka, F. Theoretical


35


Page 36 of 38

Study of the MoS2 (100) Surface: A Chemical Potential Analysis of Sulfur and Hydrogen Coverage. J. Phys. Chem. B 2000, 104, 11220–11229. (37)

Cristol, S.; Paul, J. F.; Payen, E.; Bougeard, D.; Clémendot, S.; Hutschka, F. Theoretical Study of the MoS2 (100) Surface: A Chemical Potential Analysis of Sulfur and Hydrogen Coverage. 2. Effect of the Total Pressure on Surface Stability. J. Phys. Chem. B 2002, 106, 5659–5667.

(38)

Jean-Francois Paul, E. P. Vacancy Formation on MoS2 Hydrodesulfurization Catalyst: DFT Study of the Mechanism. J. Phys. Chem. B 2003, 107, 4057–4064.

(39)

Sholl, D.; Steckel, J. A. Density Functional Theory: A Practical Introduction; John Wiley and Sons, Inc.: Hoboken, 2009.

(40)

NIST Chemistry WebBook, NIST Standard Reference Database Number 69; Linstrom, P. J., Mallard, W. G., Eds.; National Institute of Standards and Technology: Gaithersburg, MD.

(41)

Chase, M. W. NIST-JANAF thermochemical tables http://kinetics.nist.gov/janaf/ (accessed Nov 12, 2017).

(42)

CRC Handbook of Chemistry and Physics, 98th ed.; Rumble, J. R., Ed.; CRC Press: Boca Raton, 2017.

(43)

Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Atoms, Molecules, Solids, and Surfaces: Applications of the Generalized Gradient Approximation for Exchange and Correlation. Phys. Rev. B 1992, 46, 6671– 6687.


36

Page 37 of 38 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(44)

Kronberg, R.; Hakala, M.; Holmberg, N.; Laasonen, K. Hydrogen Adsorption on MoS2 Surfaces: A DFT Study on Preferential Sites and the Effect of Sulfur and Hydrogen Coverage. Phys. Chem. Chem. Phys. 2017, 19, 16231–16241.

(45)

Schweiger, H.; Raybaud, P.; Kresse, G.; Toulhoat, H. Shape and Edge Sites Modifications of MoS2 Catalytic Nanoparticles Induced by Working Conditions: A Theoretical Study. J. Catal. 2002, 207 (1), 76–87.

(46)

Steinmann, S. N.; Corminboeuf, C. A Generalized-Gradient Approximation Exchange Hole Model for Dispersion Coefficients. J. Chem. Phys. 2011, 134 (4), 44117.

(47)

Nørskov, J. K.; Bligaard, T.; Logadottir, A.; Kitchin, J. R.; Chen, J. G.; Pandelov, S.; Stimming, U. Trends in the Exchange Current for Hydrogen Evolution. J. Electrochem. Soc. 2005, 152, J23–J26.

(48)

Bligaard, T.; Nørskov, J. K.; Dahl, S.; Matthiesen, J.; Christensen, C. H.; Sehested, J. The Brønsted–Evans–Polanyi Relation and the Volcano Curve in Heterogeneous Catalysis. J. Catal. 2004, 224, 206–217.


37


Page 38 of 38

TOC Graphic


38

Comprehensive Phase Diagrams of MoS2 Edge Sites Using

Recommend Documents