Atomistic Elucidation of Sorption Processes in ... - ACS Publications

Subscriber access provided by UNIV OF DURHAM

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Atomistic Elucidation of Sorption Processes in Hydrogen Evolution Reaction on a van der Waals Heterostructure Sumit Bawari, Tharangattu N. Narayanan, and Jagannath Mondal J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b01988 • Publication Date (Web): 18 Apr 2018 Downloaded from http://pubs.acs.org on April 19, 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 25 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

Atomistic Elucidation of Sorption Processes in Hydrogen Evolution Reaction on a van der Waals Heterostructure

Sumit Bawari, T N Narayanan, and Jagannath Mondal Tata Institute of Fundamental Research-Hyderabad, Sy. No. 36/P, Gopanapally Village, Serilingampally Mandal, Hyderabad-500107, India. (Corresponding authors: [email protected] and [email protected] )

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 25

ABSTRACT:

Optimization of hydrogen evolution reaction (HER) on a catalytic surface is highly essential, and an atomistic elucidation of the key steps of HER in a given catalytic surface can provide better control for tuning productivity towards a hydrogen-based economy. However, the multi-step nature of HER poses an important challenge in this regard and this is further compounded in the recently developed atomic layers based complex heterostructured HER catalysts. A multi scale approach towards this end can provide an efficient way of dissecting the likely contribution of each step. In this work, computer simulation has been used to study the first and last steps of HER on a hexagonal Boron Nitride/Graphene (hBN/G) vertical heterostructures. Using classical Molecular Dynamics(MD) simulation combined with enhanced sampling techniques, we atomistically map the complex hBN/G heterostructure for identifying location of facile physical adsorption of protons on the heterostructure and facile desorption of evolved hydrogen from the surface. Our simulation results reveal the interface region of hBN and G surfaces or specifically the hBN-edge as free energetically favorable location for both physisorption of proton and desorption of molecular hydrogen gas. An atomistic investigation of the simulation trajectory specifically

delineates

the

role

of

zigzag

boron

terminated

edge

in

the

active

adsorption/desorption process. Subsequent density functional theory based simulations for electron transfer reaction on these classically identified adsorption/desorption sites validate their potential as facile locations for HER. Overall, our current simulation result, together with our recent experimental based findings puts forward a comprehensive multi scale picture of complete HER process on a hetero-structure and provides atomistic insights potentially useful for better material design.


2

Page 3 of 25 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. INTRODUCTION:

Hydrogen economy depends on facile hydrogen evolution reaction (HER) from substrate surface. A major thrust in this research is being expended on optimizing the substrate for catalyzing HER.1 Our recent study shows that, stacked van der Waals (vdW) vertical heterostructures (junctions) of catalytically inactive hexagonal Boron Nitride (hBN) and graphene (G), (hBN/graphene (G)) can be made active via vdW stacking.2 And these can act as efficient catalysts for HER. A similar effect has been observed from graphene/MoS2 vertical heterojunctions, where photo-electrocatalytic HER properties of the vdW heterostructure with graphene on MoS2 has the highest activity, compared to individual layers and also MoS2 on graphene.3 While these observed phenomena are fundamentally quite interesting as a complex mechanism emanating from quantum hetero-junctions of atomic layers, these are very promising from applied aspect too, as it opens possibility for designing new HER catalysts without precious metal or devoid of tiresome atomic doping processes. Though it has been found that even turbostratic stacking of G and hBN can lead to enhanced catalytic activity, the resultant structures are quite complex in atomistic dimensions, leading to a convoluted nature of HER process from these structures.2 In this regard, an atomistic understanding of the origin of HER process from these structures can lead to the designing of new catalyst materials and structures via simple stacking processes. However, as illustrated by the schematic in Figure 1, the HER process is multi-step in nature and the constituent steps are diverse both in length scales and time scales. A comprehensive understanding of both physical adsorption/desorption processes and heterogeneous charge transfer mechanism, leading to complete HER process, clearly needs adoption of a multi scale approach.4 Such a detailed approach on atomic layers based vdW structures has so far not been conducted and is highly essential for tuning their HER activity to


3


Page 4 of 25

make them comparable with benchmarked metal/metal oxide based catalysts. The current work, in combination with a recent work by Bawari et al, takes a key step towards achieving a multi scale frame-work in HER.2

FIGURE 1: Schematic for Hydrogen Evolution Reaction (HER), at different length and time scales. With adsorption and desorption limited steps between the solution and interface, at nanosecond timescales. And charge transfer limited steps between the interface and surface (S) at femto second timescales.

The mechanism of the HER is generally hypothesized to comprise of multiple steps.5 As shown in figure 1, these involve physi-adsorption of protons on the catalyst surface, electron transfer steps leading to conversion of physi-adsorbed proton to hydrogen gas molecule and finally physi-desorption of the hydrogen gas molecules from the catalyst surface. The electron


4

Page 5 of 25 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


transfer steps leading to conversion of physi-adsorbed proton (H+ or hydronium (H3O+)) to hydrogen gas molecule (H2(g)) are termed as Volmer-Heyrovsky and Volmer-Tafel mechanisms, and generally macroscopic free energetics involving these electron transfer steps are considered as an empirical metric for tuning the material towards electrocatalytic activity.1 However, it is also believed that electron transfer step notwithstanding, the steps immediately before (physiadsorption of proton on the material surface) and steps immediately after the electron transfer step (physi-desorption of hydrogen molecules from the catalytic surface) are very crucial for complete execution of the hydrogen evolution reaction.6 This is primarily because the protons need to first physically adsorb on the catalyst surface, before it can get reduced to Hydrogen gas (H2(g)) via electron transfer. Subsequently, the H2(g), obtained via electrocatalysis, needs to be physically desorbed from the surface material. This is more relevant for a complex vdW stacked hetero-structure, where it is difficult to specifically identify the role of different location(s) of the surface. In this regard, an atomistic insight on each of these steps can provide a better control towards designing better catalysts for efficient and low-cost HER. In a recent article, Bawari et al, have experimentally shown that defect free hBN/graphene heterostructures surface can catalyze the HER quite efficiently compared to its constituent defect free materials.2 In that work, Bawari et al also employed a Density functional theory (DFT) based calculation to computationally map the favorable atomistic site for facile electron transfer reaction. However, in that work, the location was mainly explored for facile electron transfer reaction by computational scanning relying on methods of trial and intuition. The current article puts together a more systematic approach of computationally simulating HER by first providing an atomistic account of the first and last steps of the hydrogen evolution reaction on hBN/graphene vdW stacked heterostructure using simulation, and subsequently computationally


5


Page 6 of 25

testing the MD-derived key sorption locations for facile electron transfer reaction. In this regard, first we use classical Molecular Dynamics (MD) simulation to computationally elucidate the physi-adsorption of H3O+ ions on the heterostructure and physi-desorption process of H2 molecule from the heterostructure for hBN/graphene surface. Since sorption steps do not require the process of bond-breaking and bond-making, Classical MD simulation acts as the suitable method of choice to simulate these physical processes using the classical approach. MD provides us with the advantage of longer time scale and larger length scale exploration in presence of explicit solvent and ions compared to DFT-based techniques. Towards this end, we employ a combination of equilibrium MD simulation and enhanced sampling method on a full classical force field and quantify the underlying free energy landscapes and associated kinetics of physiadsorption process of H3O+ ions on to the surface and physi-desorption process of H2 gas molecules from the heterostructured surface independently. This study reveals an atomistic map of preferred locations of physiadsorption of H3O+ and physidesorption of H2 gas molecules from a complex heterostructure of hBN/graphene stacked heterostructure. Subsequently, Density Functional Theory (DFT) based technique is employed on these preferred location of physiadsorption/desorption for verifying whether these are facile for efficient electron transfer reaction required for HER. Specifically, our results point towards the crucial role that interface region of hBN/graphene heterostructure can play in guiding the associated physical processes preceding and following electron-transfer process of HER and together with our previous experimental/computational venture investigating the electron-transfer process of HER on the same surface, put forward an multiscale perspective of HER.


6

Page 7 of 25 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


2. COMPUTATIONAL METHODS 2.1 Simulation Model

FIGURE 2: System considered for MD simulations, with hBN, graphene, H3O+, Cl- and H2 molecules in a water box of 5.5x5.5x4 nm3. The collective variables (CVs) used in umbrella sampling simulations [xy and z] are shown with 3 different regions on the heterostructure shown, hBN center(red), hBN edge(blue) and graphene(orange). Our system for simulation consists of a planar sheet of Graphene (4x4 nm2), stacked with a smaller hBN (1.25x1.25 nm2) in AA configuration. The choice of smaller hBN on larger


7


Page 8 of 25

graphene is to mimic our earlier experimental work on hBN/G heterostructure dimensions.7 The mutual distance of separation of hBN and graphene sheet is 0.33 nm, which is optimized before hand and also observed in experiments.8 The hBN sheet is hydrogen terminated and carbon atoms of graphene model uncharged sp2 hybridized carbon atoms and parameters are chosen from literature where the flow of water through nanotubes is studied.9 The parameters for modeling hBN and graphene sheets are given in table 1.9 A rectangular simulation box of dimension 5.5x5.5x4 nm3 is used and periodic boundary is employed in all three dimensions. The simulation box is first solvated with water molecules modeled by SPC/E classical water model.10 In the current simulation, protons are modeled by hydronium ions (H3O+). To simulate the electrochemical environment, sufficient number of H3O+ and Cl- ions pertaining to a ~0.5 M H2SO4 solution [pH = 0], are added to the solvated simulation box. One hydrogen molecule is also added to the simulation box, to study the desorption process of H2 gas molecule. The classical parameters modeling H3O+, Cl- and H2(g) are obtained from literature and tabulated in table 1.9,11 In summary, the overall system is composed of an hBN on a graphene sheet, surrounded by ~3500 water molecules and 66 H3O+ and Cl- ions each, and one H2 molecule. The non

bonding

potential

= 4

between

where = √ and ! =

−

"

any

+

two

atoms

A

and

B

is

given

as

...(1)

, between any ε , σ and q given in table 1.

TABLE 1: Force field parameters used for graphene, hBN, H3O+, Cl-, and H2O 9,11


8

Page 9 of 25 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


ε(kJ/mol)

σ(nm)

q(coulomb)

Carbon

0.236366

0.321443

0.000

center

0.3971

0.3453

0.37

edge

0.3971

0.3453

0.32

center

0.6060

0.3365

-0.37

edge

0.6060

0.3365

-0.19

H-H

0.0683

0.2813

0.00

B-H

0.0683

0.2813

-0.32

N-H

0.0683

0.2813

0.19

H3O+

0.0418

0.0801

0.4606

H2O

0.00

0.00

0.4238

H3O+

0.7732

0.3165

-0.3819

H2O

0.6501

0.3165

-0.8476

Chlorine

0.5200

0.3780

-1.00

Atom/Parameter

Boron

Nitrogen

Hydrogen

Oxygen

2.2 Simulation Method: Molecular Dynamics To investigate the free energetics and kinetics of underlying physisorption process of H3O+ ion on the hBN/graphene surface and physidesorption process of H2(g) molecules, we perform both equilibrium MD simulation and umbrella sampling simulation. While, equilibrium MD provides a kinetic picture of the overall process, umbrella sampling approach has the advantage of quantitatively computing the free energy landscape and corresponding barriers and minima related to both physisorption of H3O+ ions and physidesorption of H2 gas from the surface. Both the approaches supplement each other to atomistically map regions of facile adsorption and desorption.


9


Page 10 of 25

All MD simulations are performed using GROMACS 5.0.1 package.12 First, the electrochemical system described earlier is simulated in equilibrium condition for 1 microsecond. A time step of 2 femtosecond was used. The verlet cutoff scheme was used with a LJ cutoff of 1 nm. The electrostatic potential is treated using particle mesh Ewald approach.13 The graphene and hBN sheets are treated as “freeze-group”(as implemented in GROMACS) so as to make them immobile. Nose-Hoover thermostat is used to maintain an average temperature of 300 K.14 SETTLE algorithm is used for rigid water molecules and LINCS algorithm is used to maintain rigid tetrahedral geometry of H3O+ and linear geometry of H2 gas molecules.15,16 To quantitatively map the free energy landscape of physi-adsorption and physi-desorption of H3O+ and H2(g) respectively in the entire hBN/graphene sheet, umbrella sampling scheme is employed along two orthogonal Collective Variables (CV)(figure 2), which cover both lateral and vertical directions of the stacked material.17 The two CVs are given as the distance between the center of the sheet and the molecule along xy=%& +y (referred to as xy plane) and z directions as collective variables(figure 2). For the purpose of umbrella sampling, the Gromacs program was patched with PLUMED plugin (version 2.2.3).18 We adopt harmonic restraintsU=0.5+xy ,xy − xy- . + 0.5+/ ,0 − 0- . with Kxy = 119.503 Kcal/mol/nm2 and Kz = 119.503-780.1147 Kcal/mol. The values of Kxy and Kz were chosen such that the distributions of the corresponding reaction coordinates around the desired xy0 and z0 are Gaussian in nature and there is significant overlap among adjacent windows(figure S1). The CVs sampled by umbrella sampling ranged between 0 and 1.2 nm along xy and between 0 and 1.0 nm along z at an interval of 0.1 nm in both dimensions. The overall protocol resulted in umbrella sampling a two dimensional grid-space amounting to a total of 13 X 11 = 143 windows.


10

Page 11 of 25 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


Each of the two-dimensional windows were umbrella sampled for 500 ps, thereby amounting to a total of (143 X 0.5)= ~70 ns simulation time. The sampling schemes were repeated for both H3O+ and H2.

Finally, the underlying potential of mean force or free energy along two

dimensions was obtained using Weighted Histogram Analysis Method (WHAM).19,20 2.3 Simulation Method: Density Functional Theory Geometry relaxations and electronic structure calculations in the current article have been performed using Density Functional Theory (DFT) using the Spanish Initiative for Electronic Simulations of Thousands of Atoms (SIESTA 4.0) package.21 The system studied is a stacked hetero-structure of h-BN/Graphene (34+96 = 130 atoms). The stacking arrangement AA,is considered, similar to the MD simulations. The hetero-structure label AA corresponds to both Boron and Nitrogen atoms of h-BN (34 atoms) coinciding with Carbon atoms of the Graphene (96 atoms) sheet. As the hetero-structures remain bound due only to the van der Waals interactions between the sheets we have used pseudo-potentials generated using the norm-conserving Troullier-Martins scheme incorporating the van der Waals ‘vdW-DF2’ correction in the existing Generalized Gradient

Appoximation

(GGA)

Perdew-Burke-Ernzerhof

(PBE)

exchange

correlation

functional.22-24 The energy cutoff for the real space grid used to represent the density was set at 525 Ry for the composite heterostructure. The Broyden method was employed in a variable- cellrelaxation scheme for geometry relaxation until the maximal forces on each relaxed atom were less than 0.00001 Hartree/Bohr.25 The periodic unit cell is a cuboid with cell constant a=17.1 Angstrom with a 4x4x1 k-points sampling of the Brilliouin zone using Monkhorst-Pack scheme


11


Page 12 of 25

for electronic calculations.26 All the calculations were performed spin-unrestricted and we have used the Hirshfeld population scheme for charge density analysis.27,28 Finally we calculate ΔGH for the Hydrogenation(after the Volmer step) on the Surface(Su) in each of these systems, at different locations using the following equation.1 ΔGH = ESu-H – ESu – 0.5 EH2 - [TΔS + ZPE] ...(2) Where ZPE is the zero point energy correction which along with TΔS is a constant value of 0.24 eV added to get the correct ΔGH.


12

Page 13 of 25 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


3. RESULTS AND DISCUSSION: One of the major goals of the current work is to provide an atomistic mapping of the vdW hBN/graphene surface for facile adsorption of H3O+ ion and facile desorption of molecular hydrogen gas (H2(g)). Towards this end, free energy analysis is considered to be a quantitative approach of mapping the extent of favorability of the sorption of H3O+ ion and H2(g) across the surface and away from the surface. In this work, as detailed in the method section, we have employed umbrella sampling approach to delineate the desired free energetics of sorption process of H3O+ ion and H2(g) using two collective variables, which represent the distance of the molecules from the stacked surface in xy plane and from the z direction (figure 2). Figure 3A depicts the simulated multi-dimensional free energy landscape of sorption of H3O+ ion on the stacked surface covering a range of 0.0 to 1.2 nm in the xy direction and 0.0 to 1.0 nm in the z direction. As depicted in figure 2, xy = 0.0 nm is the location of the center top hBN sheet, xy = 0.8-0.9 nm represent the interface of hBN and graphene regions, while xy = 1.0-1.2 nm represents the pure graphene region. As can be seen in figure 3A, for xy = 0-0.7 nm, which corresponds to the pure hBN sheet, anything closer than z = 0.2 nm to the sheets is too repulsive for the molecules to approach for steric reason. We note that the thermal fluctuations in a system at 25 °C are of the order of ~0.59 kcal/mol. A closer look at figure 3A shows that at a range of xy = 0-0.7 nm, which corresponds to pure hBN surface, two free energy minima exist: at a very close molecule-material separation of z = 0.3 nm and at a slightly larger molecule-material separation of z = 0.6 nm. But these two free energy minima along xy = 0-0.7 nm are well separated by a high free energy barrier at z = 0.5 nm, potentially limiting the probability of adsorption of H3O+ ions at z = 0.3 nm. On the other hand, in the pure graphene region, corresponding to xy = 1.0-1.2 nm, figure 3A shows a minima at z = 0.1 nm (the graphene- H3O+


13


Page 14 of 25

ion contact point) and a broad minima at larger separation of z = 0.4-1.0 nm (solvated regime) mutually separated by a free energy barrier at an intermediate z value of 0.3 nm, which makes the approach of H3O+ ion unlikely on the graphene surface. However, interestingly, at around xy = 0.9 nm, which corresponds to the interface region of graphene-hBN stacked heterostructure, the adsorption of H3O+ ions from solvated regime to the material surface is barrier-less. The onedimensional free energy profiles projected along z, for a given xy, shown in figure 3C echo the similar conclusion that adsorption of H3O+ ion lacks any barrier at the edge of hBN sheet. Hence the free energy analysis overall predicts that the edge or interface region of graphene-hBN heterostructure H3O+ ions will be most facile for adsorption. Figure 3C and 3D quantify respectively the multidimensional free energy landscape and onedimensional projection along z, corresponding to the sorption process of H2(g) molecule on the hBN/graphene stacked surface. The free energetics of H2 molecule also show large free energy barrier at an intermediate z value in the pure hBN surface (xy = 0-0.7 nm). However, the barrier gradually reduces towards the hBN edge or hBN-graphene interface region, making the desorption process of H2 molecules more facile at around the edges. The relative free-energies encountered in Figure 3, are in the range of -1.5 to 3 kcal/mol with barriers of ~0.5-1.5 kcal/mol. However, as reflected in the unbiased MD–derived probability distributions of Figure 4, even these small differences in free-energy can considerably tune the adsorption/desorption dynamics of these molecules at various locations of the heterostructures. These barriers of ~0.5-1.5 kcal/mol make certain paths relatively less favorable, thereby favoring other locations for adsorption/desorption.


14

Page 15 of 25 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: Free Energy Surfaces along hBN/Graphene, for H3O+ molecule(A), and for H2 molecule(B). One dimensional projections of free energies along different xy for H3O+ molecule(C), and for H2 molecule(D)


15


Page 16 of 25

To obtain a more dynamic picture of this adsorption/desorption process, a micro-second long unbiased simulation of the whole system was performed (see methods). To compare with probabilities derived from umbrella-sampling free energetics, the probability distribution of locations visited by the molecules in a grid of xy and z are plotted, for the H2(g) and one H3O+ ion. The probability distributions obtained from independent equilibrium simulation correlates well with umbrella-sampled free energy landscapes described earlier. Consistently, a very high occupancy for H3O+ ion is seen at [xy = 0.9, z = 0.0], which corresponds to hBN edge. On the other hand, a similar but lower occupancy is seen for H2 around the same location [xy = 0.9, z = 0.0], also pertaining to the hBN edge.

FIGURE 4: Occupation probabilities along hBN/Graphene surface calculated from unbiased simulations, for H3O+ ion(left), and for H2(g)(right).


16

Page 17 of 25 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


Taken together, our simulation results provide an interesting insight that edges of the hBN/graphene vdW heterostructure can act as a location for facile adsorption of proton (H3O+) and desorption of H2(g). Graphene is completely neutral and does not interact much with ions in the solution. Stacking of the hBN sheet on graphene provides much needed polarity for catalysis. hBN itself is inactive for HER, and the top of the hBN sheet proves inefficient for adsorption of H3O+ ion and/or the desorption of H2(g)(fig 3A and B). Hence, in all likelihood, the hydrogen terminated hBN edge poses a fair possibility for easier movement of H3O+ and H2 molecules. This indicates the role of hBN edge as an active adsorber of the H3O+ ion and desorber of the H2(g), while the electron transfer steps can take place on the conducting graphene sheet.

FIGURE 5: Occupation of H3O+ and H2 molecule, along the hBN edges(A). Radial distribution functions calculated from the edge hydrogens of Nitrogen Zigzag(top), Armchair(left and right),


17


Page 18 of 25

Boron Zigzag(bottom), and center two hydrogens(Mid H) of Boron Zigzag(bottom) edges of hBN for H3O+(B) and H2(C) molecule The usefulness of sampling techniques lie in the choice of appropriate collective variables, which are dimensionally reduced reaction coordinate. While our choice of current collective variables provides a decent description of the underlying sorption behavior of small molecules on the catalytic surface, a limitation of using xy as a collective variable is the assumption that the hBN and graphene sheets are isotropic in x and y directions. While this is true for the region of graphene considered, the hBN sheet can have three different types of edges (figure 5): one armchair (left and right), boron terminated zigzag (bottom), and nitrogen terminated zigzag (top). While the umbrella sampling runs do not differentiate between these edges, the microsecondlong unbiased atomistic simulation can provide more atomistic insights into possible active sites along the edges. Figure 5 shows a scatter plot of center of mass position of H3O+ (blue) and H2 (pink), when these are very close to the hBN edge, in the range xy[0.7-1.0 nm] and z[0.0-0.2 nm]. We find that, although both molecules have fairly similar extent of occurrences in the armchair edges, H3O+ ions adsorb exclusively on the boron-terminated zigzag edges, while H2 has no occupation in this range. The high occupation of H3O+ at this location suggests that the H3O+ maximum in probability distributions corresponds to the boron-zigzag edge. On the other hand, for H2, the maxima may correspond to either the armchair or nitrogen-zigzag edges. A comparison of pair correlation function (figures 5B and C) of terminating hydrogen of all three edges with both H3O+ and H2, calculated from the unbiased simulation, provides a further quantitative evidence of edge-selectivity of these two molecules. For H3O+, as shown in Figure 5B, there appears a clear maxima in pair correlation function at 0.3 nm for boron-Zigzag edges, while it has very low correlation with other edges. On the other hand, for H2(g), as shown in


18

Page 19 of 25 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 5C, the radial distribution function shows very low distribution for Zigzag boron and nitrogen edge, while a small maximum exists for the armchair edges. The distribution increases further for H3O+, and decreases further for H2 (figure 5B and C (in pink)), if only the middle hydrogens of the Zigzag boron edges are considered for calculating the radial distribution function. This site preference, as elucidated by our classical atomistic simulation approach, tends to suggest that if the H3O+ adsorbs at the Zigzag boron edge and the electron transfer reaction takes place, its desorption of produced molecular H2 is very likely, as the probability of molecular hydrogen residing there is found to be very low. Comparison simulations were performed for a stand-alone hBN sheet in similar conditions(figure S2 and S3). hBN seems to play a greater role in the sorption behaviour, with an edge dependance. Probability distributions along the hBN surface in figure S2 show much diminished occupation for H3O+ at the hBN edge[xy = 0.9 , z = 0]. While for the H2 molecule, high occupation locations are more spread out, favouring confinement of the H2. The edge dependence seen in hBN/G is greatly diminished for H3O+ on the B zigzag edge and is almost completely lost for H2(figure S3). This proves that the enhanced sorption and selectivity effect, while steered by hBN, is not due to hBN alone, but due to the stacking of hBN on graphene. For reproducibility check, another one microsecond unbiased simulation is performed for the hBN/G system, which is independent of the first simulation. The analyses are compared for two independent trajectories in figure S4-S6. An overall agreement in the results of these independent simulations, provides more confidence in the data. Previous DFT calculations had shown that the locations on the graphene side that lie directly over a Nitrogen of Boron, yield lowest ΔGH [1.06 eV](1 eV= 23.061 kcal/mol) observed in the


19


Page 20 of 25

graphene/hBN system.2 While classical MD results show the hBN edge locations as efficient adsorbers of H3O+ ion and H2(g), which lie on the opposite side. Experiments suggest G/hBN stacks as better catalysts for HER than hBN/G stacks, but hBN/G systems also show much enhanced activity. This has also been confirmed through DFT calculations on certain locations on the hBN edge which also show medium ΔGH [~1.10 eV], that is certainly lower than pure Graphene[1.18 eV].

FIGURE 6: Occupation of H3O+ ion and H2(g) molecule, along the hBN edges highlighting Bz(Boron zigzag), Nz(Nitrogen zigzag) and Ac(Armchair) locations(left). ΔGH calculated for the same locations(shown in blue).


20

Page 21 of 25 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


Since all ΔGH encountered in the G/hBN system have been found to be positive, this implies that the formation of the S-H bond (Volmer step) is the unfavorable step.2 Locations of high H3O+ ion adsorption can help overcome this problem, but a low ΔGH is also desirable for those locations. Towards this end, in the current work, we quantum-mechanically evaluate the three locations [Boron zigzag, Nitrogen zigzag and Armchair] capable of facile sorption, for their efficiency towards electron-transfer reaction. In this regard, DFT calculations are employed and ΔGH as described in equation 2 are calculated for these locations. Our DFT-based results show that Boron zigzag and Armchair locations result in medium ΔGH [1.09 and 1.11 eV], as observed in previous studies, while Nitrogen zigzag location yields slightly higher ΔGH [1.15 eV]. Coupling with the high H3O+ ion adsorption tendency obtained via classical MD simulations, and a low ΔGH (considering the hBN/G surface) evaluated via DFT based approach, these results suggest that the Boron zigzag site is the catalytically active site for HER when hBN is stacked on Graphene. 4. CONCLUSION: Understanding HER from materials surface at an atomic level is highly desired for better catalyst design. In this regard, an atomistic mapping of the complex catalytic surface for facile execution of multi-step process can provide key insights. The current work, using classical MD simulation, takes an important step in this direction by providing atomistic picture of the first and last steps of the HER on hBN/graphene vdW stacked surface. Combining equilibrium MD and umbrella sampling approach, we have atomistically mapped the hetero structure for location of facile physiadsorption process of protons and physidesorption processes of H2(g). Our simulation suggests the interface region of hBN and graphene surface or specifically the hBN-edge as a


21


Page 22 of 25

location of choice for both facile adsorption of H3O+ and desorption of H2 (g). An atomistic investigation of the simulation trajectory specifically delineates the role of zigzag boron terminated edge in active adsorption/desorption process. Subsequent DFT based quantum calculations further validate these sites’ ability to perform facile electron transfer reaction. Our current simulation results, together with Bawari et al’s experimental findings of HER on the same hBN/graphene surface, provides a complete account of essentially all steps of HER process and provides key insights on the facile location of electron transfer process.2 We believe our multi scale approach in understanding HER from atomistic perspective will open

new

opportunities for materials design. We note that, the hBN considered in the present work is Hydrogen terminated across all edges. While this is required for charge neutralization in MD simulations, the Nitrogen-terminated edges may oxidise under experimental HER conditions. We plan to probe this in future work, by an in-depth DFT study. SUPPORTING INFORMATION The Supporting Information is available free of charge on the ACS Publications website at DOI: Distribution of Umbrella sampled CV, probability distribution and RDF for stand-alone hBN, and comparison of probability counts, adsorption sites and RDF between two independent unbiased simulations.(PDF). ACKNOWLEDGMENT Authors thank the financial support received from Tata Institute of Fundamental Research, Hyderabad. REFERENCES (1)

Seh, Z. W.; Kibsgaard, J.; Dickens, C. F.; Chorkendorff, I.; Nørskov, J. K.; Jaramillo, T. F. Combining Theory and Experiment in Electrocatalysis: Insights into Materials Design. Science (80-. ). 2017, 355 (6321), eaad4998.

(2)

Bawari, S.; Kaley, N. M.; Pal, S.; Vineesh, T. V.; Ghosh, S.; Mondal, J.; Narayanan, T. N. On the Hydrogen Evolution Reaction Activity of Graphene-hBN van der Waals Heterostructures Phys. Chem. Chem. Phys., 2018, DOI: 10.1039/C8CP01020J


22

Page 23 of 25 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


(3)

Biroju, R. K.; Das, D.; Sharma, R.; Pal, S.; Mawlong, L. P. L.; Bhorkar, K.; Giri, P. K.; Singh, A. K.; Narayanan, T. N. Hydrogen Evolution Reaction Activity of Graphene–MoS2 van Der Waals Heterostructures. ACS Energy Lett. 2017, 2 (6), 1355–1361.

(4)

Skúlason, E.; Tripkovic, V.; Björketun, M. E.; Gudmundsdóttir, S.; Karlberg, G.; Rossmeisl, J.; Bligaard, T.; Jonsson, H.; Nørskov, J. K. Modeling the Electrochemical Hydrogen Oxidation and Evolution Reactions on the Basis of Density Functional Theory Calculations. J. Phys. Chem. C 2010, 114 (42), 18182–18197.

(5)

Feng, J.-X.; Wu, J.-Q.; Tong, Y.; Li, G.-R. Efficient Hydrogen Evolution on Cu Nanodots-Decorated Ni3S2 Nanotubes by Optimizing Atomic Hydrogen Adsorption and Desorption. J. Am. Chem. Soc. 2017, 140(2), 610−617.

(6)

Liang, Z.; Ahn, H. S.; Bard, A. J. A Study of the Mechanism of the Hydrogen Evolution Reaction on Nickel by Surface Interrogation Scanning Electrochemical Microscopy. J. Am. Chem. Soc. 2017, 139 (13), 4854–4858.

(7)

Vinod, S.; Tiwary, C. S.; da Silva Autreto, P. A.; Taha-Tijerina, J.; Ozden, S.; Chipara, A. C.; Vajtai, R.; Galvao, D. S.; Narayanan, T. N.; Ajayan, P. M. Low-Density ThreeDimensional Foam Using Self-Reinforced Hybrid Two-Dimensional Atomic Layers. Nat. Commun. 2014, 5, 4541.

(8)

Gao, G.; Gao, W.; Cannuccia, E.; Taha-Tijerina, J.; Balicas, L.; Mathkar, A.; Narayanan, T. N.; Liu, Z.; Gupta, B. K.; Peng, J.; et al. Artificially Stacked Atomic Layers: Toward New van Der Waals Solids. Nano Lett. 2012, 12 (7), 3518–3525.

(9)

Won, C. Y.; Aluru, N. R. Water Permeation through a Subnanometer Boron Nitride Nanotube. J. Am. Chem. Soc. 2007, 129 (10), 2748–2749.

(10)

Jorgensen, W. L.; Madura, J. D. Solvation and Conformation of Methanol in Water. J. Am. Chem. Soc. 1983, 105 (6), 1407–1413.

(11)

Jang, S. S.; Molinero, V.; Çaǧin, T.; Goddard, W. A. Effect of Monomeric Sequence on Nanophase-Segregated Structure and Water Transport in Nafion 117. ACS Div. Fuel Chem. Prepr. 2003, 48 (2), 510–511.

(12)

Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindah, E. Gromacs: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1–2, 19–25.

(13)

Sagui, C.; Darden, T. A. MOLECULAR DYNAMICS SIMULATIONS OF BIOMOLECULES: Long-Range Electrostatic Effects. Annu. Rev. Biophys. Biomol. Struct. 1999, 28 (1), 155–179.


23


Page 24 of 25

(14)

Evans, D. J.; Holian, B. L. The Nose–Hoover Thermostat. J. Chem. Phys. 1985, 83 (8), 4069–4074.

(15)

Miyamoto, S.; Kollman, P. A. Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13 (8), 952–962.

(16)

Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. {LINCS}: A Linear Constraint Solver for Molecular Simulations. J. Comp. Chem. 1997, 18, 1463–1472.

(17)

Torrie, G. M.; Valleau, J. P. Nonphysical Sampling Distributions in Monte Carlo FreeEnergy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23 (2), 187–199.

(18)

Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New Feathers for an Old Bird. Comput. Phys. Commun. 2014, 185 (2), 604–613.

(19) Grossfield, Alan, "WHAM: the weighted histogram analysis method", version 2.2.3, http://membrane.urmc.rochester.edu/content/wham (20) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I: The Method. J. Comput. Chem. 1992, 13(8), 1011−1021. (21) Soler, J. M.; Artacho, E.; Gale, J. D.; García, A.; Junquera, J.; Ordejón, P.; SánchezPortal, D. The SIESTA method for ab initio order-N materials simulation. J. Phys. Condens. Matter 2002, 14(11), 2745. (22) Troullier, N.; Martins, J. L. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B 1991, 43(3), 1993. (23) Lee, K.; Murray, É. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Higher-accuracy van der Waals density functional. Phys. Rev. B 2010, 82(8), 081101. (24) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77(18), 3865. (25) Johnson, D. D. Modified Broyden’s method for accelerating convergence in selfconsistent calculations. Phys. Rev. B 1988, 38(18), 12807. (26) Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13(12), 5188. (27) Hirshfeld, F. L. Bonded-atom fragments for describing molecular charge densities. Theor. Chim. Acta 1977, 44 (2), 129–138. (28) Fonseca Guerra, C.; Handgraaf, J. W.; Baerends, E. J.; Bickelhaupt, F. M. Voronoi deformation density (VDD) charges: Assessment of the Mulliken, Bader, Hirshfeld, Weinhold, and VDD methods for charge analysis. J. Comput. Chem. 2004, 25(2), 189-210.


24

Page 25 of 25 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


Table of Contents(TOC) Graphic:


25

Atomistic Elucidation of Sorption Processes in ... - ACS Publications

Recommend Documents