Reversibility Iteration Method for Understanding Reaction Networks

Sep 9, 2016 - School of Chemistry and Chemical Engineering, The Queen's University of ... In this paper, a new approach, the so-called reversibility i...
0 downloads 0 Views 2MB Size
Subscriber access provided by Northern Illinois University

Article

Reversibility Iteration Method to Understand Reaction Networks and to Solve Micro-Kinetics in Heterogeneous Catalysis Jian-Fu Chen, Yu Mao, Hai-Feng Wang, and Peijun Hu ACS Catal., Just Accepted Manuscript • DOI: 10.1021/acscatal.6b02405 • Publication Date (Web): 09 Sep 2016 Downloaded from http://pubs.acs.org on September 9, 2016

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 free 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 accessible to all readers and 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.

ACS Catalysis 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 31

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

ACS Catalysis

Reversibility Iteration Method to Understand Reaction Networks and to Solve Micro-Kinetics in Heterogeneous Catalysis Jian-Fu Chen1, Yu Mao2, Hai-Feng Wang1*, and P. Hu1,2* 1

Key Laboratory for Advanced Materials, Research Institute of Industrial Catalysis and Centre for Computational Chemistry, East China University of Science and Technology, 130 Meilong Road, Shanghai, P.R. China, 200237. Email: [email protected]

2

School of Chemistry and Chemical Engineering, The Queen’s University of Belfast, Belfast, BT9 5AG, UK. Email: [email protected]

Abstract Solving micro-kinetics of catalytic systems, which bridges microscopic processes and macroscopic reaction rates, is currently vital for understanding catalysis in silico. However, traditional micro-kinetic solvers possess several drawbacks that make it slow and unreliable for complicated catalytic systems. In this paper, a new approach, the so-called reversibility iteration method (RIM), is developed to solve micro-kinetics for catalytic systems. Using the chemical potential notation we previously proposed to simplify the kinetic framework, the catalytic systems can be analytically illustrated to be logically equivalent to the electric circuit, and the reaction rate and coverage can be calculated by updating the values of reversibilities. Compared to the traditional modified Newton iteration method (NIM), our method is not sensitive to the initial guess of the solution and typically undergoes fewer iteration steps. Moreover, the method does not require arbitrary-precision arithmetic and has a 1

ACS Paragon Plus Environment

ACS Catalysis

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

higher probability to solve the system successfully. These features make it about 1000 times faster than modified Newton iteration method for the systems we tested. Moreover, the derived concept and the mathematical framework presented in this work may provide new insight into catalytic reaction networks. Keywords: micro-kinetics; heterogeneous catalysis; reversibility iteration method; turnover frequency; catalytic mechanism; catalysis design

2

ACS Paragon Plus Environment

Page 2 of 31

Page 3 of 31

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

ACS Catalysis

1. Introduction Rational design of catalysts assisted by first principles calculations is currently becoming an increasingly important topic in heterogeneous catalysis since the traditional trial-and-error method is difficult to meet the demands of the rapidly developing catalyst industry.1-12 In the design process, solving micro-kinetics of catalytic systems is one of the crucial steps; with the energetics of reaction networks calculated by density functional theory (DFT), micro-kinetics can predict the reaction rate (turnover frequency, TOF) and coverage of surface species, which are vital for determining the reactivity and selectivity for a given catalyst.13-21 According to the steady-state approximation,13-15 micro-kinetics involves the simultaneous solution of a set of nonlinear differential equations, which includes the mass balances of species under reaction conditions (temperature, pressure) and invariance of surface species coverages with respect to time:13,22,23



i

1

(1)

i

 i 0 t

(2)

where the derivative of coverage (∂θi/∂t) is determined by the rate equations of elementary reactions in the reaction network. One of the simplest and most widely used methods to solve such a set of equations is the modified Newton iteration method (NIM). Generally, the method exhibits a good convergence behavior for normal equations. Some problems, however, arise when solving micro-kinetics for practical catalytic systems that involve tens even hundreds of reactions.22-24 Firstly, 3

ACS Paragon Plus Environment

ACS Catalysis

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 4 of 31

the differential equations (2) are very “stiff” in most cases of practical catalytic systems, which means that the rate constants typically vary many orders of magnitude, and the coverage ranges of different surface species also span enormously from extremely close to zero to unity. Secondly, the initial guess of adsorbate coverages, which is usually required in the root-finding procedures such as NIM, influences considerably the solutions; unreasonable initial guesses may lead to a very slow convergence. These features may result in oscillations and erratic solutions; therefore, special numerical techniques are required. Since the standard-precision arithmetic (typically 16 decimals) is far from sufficient to obtain a stable solution, one idea is to employ

arbitrary-precision

arithmetic.22,24

However,

it

would

be

rather

time-consuming and memory-demanding when the precision is too high. Some studies have been devoted into more stable and accurate solution techniques of micro-kinetics. Gusmão and Christopher23 proposed a so-called “general and robust” quasi-Newton approach to solve micro-kinetics for catalytic systems. They found that in micro-kinetics, adsorbate mass balances can be automatically deduced by the stoichiometric matrix, and then the Taylor series expansion of the adsorbate mass balances can be performed to obtain an analytical Jacobian and Hessians matrixes that encompass first-order and second-order partial derivatives of the mass balance equations, respectively. These features ensure a good convergence and high accuracy. Kozuch et al.25,26 proposed the energetic span model that tried to estimate the TOF of a complex reaction network by the energy difference between the TOF-determining transition state and the TOF-determining intermediate, which is quite reasonable for 4

ACS Paragon Plus Environment

Page 5 of 31

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

ACS Catalysis

some organometallic reactions shown in their work. In this contribution, an effective “reversibility iteration method” (RIM) is developed to understand the reaction mechanism and to solve micro-kinetics for the catalytic reactions. The chemical potential notation27 that proposed in our previous work is employed to formulate micro-kinetic equations, and an analogy of catalytic systems to the electric circuit can be clearly illustrated by these equations. The RIM can be derived for all kinds of catalytic networks, including series reactions, parallel reactions, and complex reactions that involves both of them. Instead of solving differential equations, the reaction rate and coverage can be calculated by updating the values of reversibilities in our method. Our approach is successfully applied to two practical catalytic systems: ammonia synthesis and iodine reduction reaction (IRR). The validity, robustness, and effectiveness of our method are evaluated and its underlying reasons are discussed; the results show that our method is not sensitive to the initial guess, and about 1000 times faster than NIM.

2. Methods and construction of the iteration cycle All the rate equations in this paper are formulated within the framework of chemical potential that is comprehensively described in the supporting information (SI-1-0). For a generalized two-step heterogeneous catalytic reaction: 1.adsorption 2.desorption R(g)  n*  nI*  P(g)  n*

(3)

where R, I, and P represent reactants, intermediates, and products, respectively; n is the number of active free sites involved in the reaction (typically n would be one or 5

ACS Paragon Plus Environment

ACS Catalysis

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 6 of 31

two in common catalytic reactions). Under the steady-state condition, the rate equation of reaction (3) can be eventually formulated as: R

k T e RT (1  zt ) n r  B   ,o  2 ,o * 1 h e RT  e RT

(4)

where:

* 

1 1  n z1 K eq1

(5)

P  R

zt  e

z1 

e

RT

 z1 z2

 2 ,o

1 ,o

RT

 zt e RT

1 ,o

 2 ,o

(6)

(7)

e RT  e RT  R   Io

K eq1  e

RT

(8)

Detailed derivations of above equations can be found in SI-1-2. z1, z2, and zt are the reversibilities28,29 of the adsorption, desorption, and total reaction, respectively; θ* represents the coverage of active free sites; µR and µP stand for the chemical potential of reactants and products, respectively, while µIo, µ1≠,o, and µ2≠,o are the standard chemical potentials of the intermediate, the transition state of the adsorption process, and transition state of the desorption process, respectively. Then, rate equation (4) can be rearranged as:

6

ACS Paragon Plus Environment

Page 7 of 31

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

ACS Catalysis

r

=

k BT h

1  zt e

1 ,o - R

 2 ,o - R

RT

RT

e

*n =

1  zt h e k BT

1 ,o - R RT

*

-n

h  e k BT

 2 ,o - R RT

 *- n

(9)

U RR  RP

where:

U 1zt

(10)

U can be defined as the driving force of the reaction, which is determined by the reversibility. When the reaction is far from the equilibrium, the reversibility will be small, corresponding to a large value of U. RR and RP in equation (9) are composed of the kinetic terms (RK) and thermodynamic terms (RT):

RK ,R 

1 1 ; RK ,P  kR kP

(11)

RT  * n  (1  n z1Keq1 )n

(12)

RR  RK,R RT ; RP  RK,P RT ;

(13)

The k parameter in the kinetic term, equation (11), is the rate constant of the corresponding elementary reaction referring to gaseous reactant R, defined as:  ,o

 ,o

k T - GR k T - GP kR  B e RT ; kP  B e RT h h

(14)

GR  1,o  R ; GP  2,o  R

(15)

Substituting equations (10)-(15) into equation (9), we can obtain the simplified equation r = U/(RR + RP). Interestingly, this equation is exactly the same form as that of Ohm's law in electronics: The driving force of the reaction, U, can be considered as 7

ACS Paragon Plus Environment

ACS Catalysis

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 8 of 31

the “reaction voltage” and RR and RP represent the “resistances” that include “kinetic resistance” and “thermodynamic resistance”, while r is the “reaction current”. Figure 1a shows the energy profile and the corresponding electric circuit diagram of the reaction. After carefully examining this formulation for gas-phase reactions and more complex networks involving both series and parallel reactions (Fig. 1b) in heterogeneous catalysis (SI-1-1 to SI-1-5), we find that all catalytic systems can be expressed in the above manner, namely in terms of the kinetic resistance, thermodynamic resistance, and driving force. Specifically, we construct a complex heterogeneous reaction (Fig. 1b) to delineate the results. It is a typical heterogeneous reaction that contains multiple reaction pathways with shared and divided pathways, in which rate equations can be eventually formulated as: R kBT RT r e (1  zt )*n ( h

1  I0,o

e

RT

)

(16)

1



1 e

 R1,o

 P1,o

RT

RT

e

1

 e

 R ,2o

 P ,2o

RT

 e RT

It can be simplified by the “resistance term” we defined before:

r

U RI 0  RI

1 1 1   RI RR1  RP2 RR2  RP2

(17)

(18)

It is noticeable that the total resistance here is in accordance with the way one calculates the electric resistance in an electric circuit. That is to say, the total resistance of series elementary reactions is the sum of their resistances, while the 8

ACS Paragon Plus Environment

Page 9 of 31

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

ACS Catalysis

reciprocal of the total resistance of parallel elementary reactions can be calculated by adding the reciprocals of the resistances of each elementary step. To show the applicability of our method in practical conditions, we chose hydrogen evolution reaction (HER), involving both series and parallel reactions, as an example (see details in SI-2-1). More generally, for a serial catalytic reaction (Fig. 1c):

R0  n*  nR1*  nR2* ...  nRm1*  Rm  n *

(19)

Defining the kinetic and thermodynamic resistances as we did before:

 ,o i

 i

G  

RK ,i 

kT   R ; ki  B e h

Gi ,o RT

1 1 ; RT  n ; Ri  RK ,i RT ki *

(20)

(21)

then the total resistance Rtot is the sum of individual resistances, and the overall driving force of the reaction (voltage) is:

U  z0  zt

(22)

where z0 = 1. Specifically, the voltage and reversibility on each resistance can be deduced as: j  i 1

Ui 

j i

 z  z j

j 0

j 0

j i

zi 

j



RK ,i RK ,tot

U

(23)

jm

zt  j 1 RK , j   j i 1 RK , j j  i 1

j m

zt  j 1 RK , j   j i RK , j

The final reaction rate is (the details of above derivations can be found in SI-1-6):

9

ACS Paragon Plus Environment

(24)

ACS Catalysis

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 10 of 31

R

k T e RT (1  zt ) n U U r B * = *n = i h RK ,tot Rtot m i1 e RT

(25)

where: 0  io j i

i  m 1 1 *

  1



n

e

RT

i 1

z

(26)

j

j 1

The final reaction rate is (the detailed derivation can be found in SI-1-6): R

k T e RT (1  zt ) n U U n r B  =  =  * * i h RK ,tot Rtot m RT e i1

(27)

The formulation of electric circuit-like rate equations gives rise to an interesting understanding on catalytic processes. It should be noted that similar concepts and equations have appeared in the literature.26,30-33 We expands the concepts further based on a much more complicated reaction system (with multiple reaction sites (n*)) and different kinetic framework (with reversibility and free site). The generalized model ( R n (g)  n*  nI*  Pn (g)  n * )

captures

the

core

character

of

many

homogeneous/heterogeneous catalytic reactions, which includes gas phase reactions (n=0), simple heterogeneous27 or homogeneous catalytic reactions (n=1), and complex heterogeneous catalytic reactions (n≥2) with dissociation and association of polyatomic molecules, which depicts the real situation better and takes multiple reaction sites (n*) into consideration, resulting in a far more complex rate equations than simple linear algebra. In this work, we derived the circuit-like catalytic diagram in a new approach. Although one may explain kinetic results using different kinetic

10

ACS Paragon Plus Environment

Page 11 of 31

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

ACS Catalysis

frameworks, our definition of the resistance in equation (11)-(15) provides a clear physic picture of the origin of resistances, which is a result from high barrier (kinetic resistance) when the adsorption is too weak, or low free site (thermodynamic resistance) when the adsorption is too strong (see further discussions below). More importantly, we developed the reversibility iteration method (RIM) based on above formulation to efficiently solve micro-kinetics by iterating reversibility and coverage terms. The flow chart of RIM is illustrated in Figure 2. Following the steps on Figure 2, we shall firstly formulate and calculate the values of kinetic resistances using equations (20) and (21) after constructing all the elementary reactions (step 1). Then, the rate equations should be changed to the form of “standard serial reaction” to obtain the revised expression of coverage and equivalent kinetic resistance (R’K,i) with reversibilities (zi), which are dependent on reaction mechanisms (step 2). Assuming the initial values of reversibilities to be unities or random values (step 3), we can calculate the values of coverage and R’K,i (step 4), and update the values of the reversibilities by substituting the values of R’K,i into equation (24) (step 5). Iterating from step 3 to 5 until the difference of reversibilities reaches the given tolerance (step 6), we can obtain the coverages and reaction rates with the obtained zi and Rk (step 7). We will show two applications of our approaches for solving the micro-kinetics of practical systems, ammonia synthesis and iodine reduction reaction (IRR), to demonstrate that our approach can replicate the results of two previously reported systems.

3. Results and discussion 11

ACS Paragon Plus Environment

ACS Catalysis

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 12 of 31

2NH3 ) is a very important catalytic reaction in Ammonia synthesis ( N2  3H2 

industry, the mechanism of which has been carefully examined34-37 and can be expressed as follows: H 2  2 *   2H * N 2  2 *   2N * N *  H *   NH *  * NH *  H *   NH 2 *  * NH 2 *  H *   NH 3  2 *

Following the steps in Figure 2, we can firstly write all the rate equations (step 1) in terms of the rate constants, coverages of intermediate, and reversibilities using the chemical potential notation:

kT r1  B e h kT r2  B e h

kT r3  B e h kT r4  B e h kT r5  B e h

H2  1 ,o

*2 (1  z1 )

(28)

*2 (1  z2 )

(29)

 N H (1  z3 )

(30)

NHH (1  z4 )

(31)

 NH  H (1  z5 )

(32)

RT

 N2  2 ,o RT

N  H  3 ,o RT

NH  H 4 ,o RT

NH2  H  5 ,o RT

2

The surface coverages of these intermediate species, θi, can be expressed in terms of zi according to the definition of reversibility:  H 2  2 H *

H  e

2 RT

z1 *  M 1 z1 * 12

ACS Paragon Plus Environment

(33)

Page 13 of 31

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

ACS Catalysis

 N2  2 N *

N  e

 NH  e

 NH  e

z2* 

2 RT

1 1  N  H   NH * 2 2 2 2 RT

1  N  H2   NH 2* 2 2 RT

2

z2*

M2

(34)

z1 z2 z3*  M 3 z1 z2 z3*

(35)

z1 z2 z3 z4*  M 4 z1 z2 z3 z4*

(36)

where Mi is the corresponding equilibrium constants referring to the gaseous reactants. Since  H   N   NH   NH 2 +* =1 , we have:

*1  1  M1 z1  M 2 z2  M 3 z1 z2 z3  M 4 z1 z2 z3 z4

(37)

Then, the rate equations can be simplified and rearranged as:

kT r1  B e h kT r2  B e h 1

k T 2 r3  B e h 1

k T 2 r4  B e h

r5 

k BT e h

H2  1 ,o RT

*2 (1  z1 ) 

N2  2 ,o

1 2 RT

RT

*2 (1  z2 ) 

*2 (1  z1 )

*2 (1  z2 )

1 3  N   H   5 ,o 2 2 2 2 RT

 *2 z1 z 2 (1  z3 )

 *2 z1 z 2 (1  z3 ) 

 N 2   H 2   4 ,o

 *2 z1 z 2 z3 (1  z 4 ) 

(39)

RK ,2

 N 2   H 2   3 ,o

RT

(38)

RK ,1

RK ,3

 *2 z1 z 2 z3 (1  z 4 )

 *2 z13 z 2 z3 z 4 (1  z5 ) 

RK ,4

 *2 z13 z 2 z3 z 4 (1  z5 ) RK ,5

(40) (41)

(42)

where equations (20) and (21) are used to convert the pre-exponential and chemical potential terms into the kinetic resistance terms. According to the mass balance law, we can obtain the following relations on the reaction rates and reversibilities:

13

ACS Paragon Plus Environment

ACS Catalysis

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

r

2 r1  2 r2  r3  r4  r5 3

zt 

z13 z 2 z3 z 4 z5

Page 14 of 31

(43) (44)

Aiming to employ the iteration method we previously derived from the general serial catalytic reaction, we can modify the rate equations to balance the stoichiometric numbers: 3 *2 (1  Z1 ) 2 1 2 1  z1 2 3 r  r1   (1  z )  * 1 3 RK ,1 3 (1  z1  z12 ) RK' ,1 ' 1

r2'  2r2 

(45)

 2 Z (1  Z 2 ) 1 (1  z2 ) 2 3 2 * z1 (1  z2 )  * 1 ' RK ,2 RK ,2 z13

(46)

 2 Z Z (1  Z 3 ) 1 1 2 3 * z1 z2 (1  z3 )  * 1 2' RK ,3 z1 RK ,3

(47)

r3'  r3 

*2 Z1Z 2 Z 3 (1  Z 4 ) 1 1 2 3 r  r4  * z1 z2 z3 (1  z4 )  RK ,4 z1 RK' ,4

(48)

 2 Z Z Z Z (1  Z 5 ) 1 2 3 * z1 z2 z3 z4 (1  z5 )  * 1 2 3' 4 RK ,5 RK ,5

(49)

' 4

r5'  r5 

where Z1 = z13 ; Z 2 = z 2 ; Z 3  z3 ; Z 4  z 4 ; Z 5  z5

(50)

3

' K ,1

R

z1 3 (1  z1  z12 ) '  RK ,1 ; RK ,2 =RK ,2 ; RK' ,3  RK ,3 z1; RK' ,4 =RK ,4 z1 ; RK' ,5  RK ,5 3 2 1  z1 2(1  z2 )

(51)

With the updated Rk,i’ and Zi, we completed step 2 in Figure 2: the rate equations of ammonia synthesis are rearranged into the form of “standard serial reaction”. Next, initializing all Zi to be unity (step 3) and calculating the chemical potentials with

14

ACS Paragon Plus Environment

Page 15 of 31

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

ACS Catalysis

quantum-chemical calculations, the values of Rk,i’ can be determined by substituting them into equation (51) (step 4). Then, the values of Zi, as stated in step 5, can be updated with the values of Rk,i’ using equation (24), thus forming an iteration cycle. Eventually, we will obtain Zi iteratively when the variation of their values between two iterations is smaller than the tolerance, which can also be deemed as a self-consistent method. After obtaining the values of zi, the free site ( * ) can be determined using equation (36), and the final reaction rate/turnover frequency (TOF) is:

r

U *2  R  K,i

(52)

where U 1 zt , the driving force of the reaction. Similar procedures are also employed to solve iodine reduction reaction (IRR), which is a very important reaction in dye-sensitized solar cells (DSSC)5,38-40 (see SI-2-2 for details). Then, we thoroughly examined the reliability and behavior of the RIM and made a comparison of our approach with NIM, the results of which will be displayed below. It should be noted that the micro-kinetic data here for ammonia synthesis and iodine reduction reaction are taken from the work of Wang et al.41 and Wang et al.38, respectively. Also, it should be mentioned that an automatic implementation code based on MATLAB is used to obtain the steady-state solution as the benchmarks (used to evaluate the relative error) by modified Newton’s method with various initial guess strategies and arbitrary precision arithmetic (see SI-2 for details). For the ammonia synthesis, the variation of the TOF along with the adsorption energy 15

ACS Paragon Plus Environment

ACS Catalysis

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

o

Page 16 of 31

o

of N (μN ), i.e. the volcano curve, was obtained by setting μN as a variable (Fig. 3a). o

We can see that the optimal μN is about -0.6 eV, which agrees well with the previous work in the literature.34 Aiming to verify the reliability of RIM, we calculated the relative error of turnover frequency (TOF), △r/rb (subscript b indicates benchmark; o

△r=|r- rb|, where r is the calculated TOF with either NIM or RIM), for each μN with the tolerance setting being 1×10-14, the results of which are displayed in Figure 3b. It is clear that most of relative errors are under this level. Moreover, all the errors are smaller than 1×10-13, which ensures the reliability of our method. Figure 3c displays the number of iterations on different adsorption energies; interestingly, on the two sides of the volcano curve, our approach converges very fast while the number of iterations soars near the peak of the volcano curve. This is reasonable because on the two sides of the volcano curve, either adsorption or desorption process will dominate the reaction network (rate-limiting) and make the differential equations easier to solve. The results for the iodine reduction reaction are illustrated on Figure 4. In Figures 4a and 4b, we set two tolerances (1×10-10 and 1×10-5) to make a comparison, illustrating that all the relative errors are smaller than the tolerance. In addition, Figures 4c and 4d give rise to the same trend for the number of iterations on different adsorption energies as that in Figure 3c; the method usually converges very fast except on a small o

range of μI near the peak of the volcano curve. The results above demonstrate that the RIM is effective. In order to further access its efficiency, a comparison between the performance of the RIM and that of the NIM for solving the micro-kinetics was carried and the result is illustrated in Figure 5. As the 16

ACS Paragon Plus Environment

Page 17 of 31

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

ACS Catalysis

iteration methods always need and rely on initial guesses, we generate 100 random initial guesses to test the performance of RIM and NIM at different adsorption energies along the volcano curve. As shown in Figure 5a, in which the blue line indicates the probability (P) of obtaining converged and correct results from RIM and the black and red lines represent the probability of getting (i) converged but wrong results; and (ii) converged and correct results from NIM, respectively. We find that at every point, at least in this case, the P for blue line is always 100%. The performance of NIM is relatively unstable and not reliable for some certain adsorption energies. Therefore, we arbitrarily chose two adsorption energies (0 and 2 eV) to investigate its influence on the solution. We can see from Figure 5b that for the adsorption of 2 eV (black line), the probability indeed increases as a higher precision was set, whereas for 0 eV (red line) the probability hardly changes. Therefore, setting a higher precision cannot ensure a successful solution from NIM, while RIM is always reliable in this case (for both 0 eV and 2 eV). From Figure 5c, we can see that the convergence speed of RIM to the “benchmarks” is much faster than NIM and gives rise to a reasonable solution even within a few steps (