Chinese Journal of Chemical Physics  2018, Vol. 31 Issue (4): 485-491

#### The article information

Yangyunli Sun, Shuo Zhang, Wen-hua Zhang, Zhen-yu Li

Theoretical Study of Adsorption and Dehydrogenation of C2H4 on Cu(410)

Chinese Journal of Chemical Physics, 2018, 31(4): 485-491

http://dx.doi.org/10.1063/1674-0068/31/cjcp1805120

### Article history

Accepted on: June 12, 2018
Theoretical Study of Adsorption and Dehydrogenation of C2H4 on Cu(410)
Yangyunli Suna, Shuo Zhanga, Wen-hua Zhangb,c,d, Zhen-yu Lia,c
Dated: Received on May 28, 2018; Accepted on June 12, 2018
a. Hefei National Laboratory of Physical Sciences at the Microscale, University of Science and Technology of China, Hefei 230026, China;
b. CAS Key Laboratory of Materials for Energy Conversion, Department of Materials Science and Engineering, University of Science and Technology of China, Hefei 230026, China;
c. Synergetic Innovation Center of Quantum Information Quantum Physics, University of Science and Technology of China, Hefei 230026, China;
d. Department of Applied Mathematics, School of Physics and Engineering, Australian National University, Canberra, ACT 2600, Australia
*Author to whom correspondence should be addressed. Wen-hua Zhang, E-mail:whhzhang@ustc.edu.cn; Zhen-yu Li, E-mail:zyli@ustc.edu.cn
Abstract: Adsorption and dehydrogenation of ethylene on Cu(410) surface are investigated with firstprinciples calculations and micro-kinetics analysis. Ethylene dehydrogenation is found to start from the most stable π-bonded state instead of the previously proposed di-σ-bonded state. Our vibrational frequencies calculations verify the π-bonded adsorption at step sites at low coverage and low surface temperature and di-σ-bonded ethylene on C-C dimer (C2H4-CC) is proposed to be the species contributing to the vibrational peaks experimentally observed at high coverage at 193 K. The presence of C2H4-CC indicates that the dehydrogenation of ethylene on Cu(410) can proceed at temperature as low as 193 K.
Key words: Dehydrogenation    Catalysis    Surface reactioin
Ⅰ. INTRODUCTION

Adsorption and reaction of hydrocarbons on transition metal surfaces are an interesting topic due to their relevance in heterogeneous catalysis and nanotechnology. Interaction between hydrocarbon and the metal substrate has effects on various processes like hydrogenation [1, 2], dissociation [3], and polymerization [4]. Recently, copper is widely used as a catalyst to grow two-dimensional graphene sheets via chemical vapor deposition of methane, ethane, or other hydrocarbon precursors [7]. Generally, copper has a much lower activity than that of main-group transition metals such as nickel [5, 6], but it is higher than silver and gold. Understanding the interaction between different copper surfaces and hydrocarbons will not only enrich the knowledge of surface science but also help to understand the mechanism of graphene growth [8, 9].

The interactions of ethylene with different copper surfaces have been investigated with a long history [10, 11]. Various experimental methods, such as X-ray absorption (XAS) and emission (XES) spectroscopies, temperature programmed desorption (TPD), high resolution electron energy loss spectra (HREELS), and infrared reflection-adsorption spectroscopy (IRAS), are adopted to characterize the adsorption configuration, electronic structure, and chemical reactivity of ethylene on Cu(111), Cu(100), Cu(110), and Cu(210) surfaces [10-13]. On these surfaces, the $\pi$-bonded ethylene adsorption configuration at top site is identified at low temperature.

Recently, the adsorption of ethylene on Cu(410) has been intensively investigated [14-16]. At low temperature and low coverage, ethylene binds with the Cu(410) surface with a $\pi$-bonded configuration. At the same time, a di-$\sigma$-bonded adsorption configuration was suggested based on HREELS spectra at higher surface temperature and/or higher surface coverage [14-16]. The di-$\sigma$-bonded ethylene was also proposed to be the precursor of the dehydrogenation reaction [15, 16]. XPS and HREELS experiments suggested that the reactive site of dehydrogenation is located at step edges [15]. The di-$\sigma$-bonded ethylene based dehydrogenation picture, however, was challenged by a recent study combining experimental and computational characterizations [17]. Density functional theory (DFT) calculations suggested that the di-$\sigma$-bonded configuration is a meta-stable state on Cu(410). A different assignment to the experimentally observed TPD peaks was also proposed [17]. Since controversies still exist for ethylene adsorption and dehydrogenation on Cu(410), a systematic study on this topic is desirable.

In this work, various adsorption configurations of ethylene on Cu(410) are investigated and dehydrogenation processes via different paths are explored. Microkinetic analysis suggests that the most possible dehydrogenation path of ethylene proceeds via ${\rm CH}_{\rm 2}$${\rm CH}_{\rm 2}$$\rightarrow$${\rm CH}_{\rm 2}CH\rightarrowCHCH\rightarrowCCH\rightarrowCC. On Cu(410), ethylene prefers to dehydrogenate directly from the most stable \pi-bonded configuration rather than the meta-stable long-bridge state. By comparing the calculated vibrational frequencies with experimental results [15], we further confirm the appearance of \pi-bonded ethylene at low coverage and low surface temperature. The experimentally observed HREELS peak at 193 K is assigned to di-\sigma-bonded ethylene on complete dehydrogenated product C-C. Ⅱ. COMPUTATIONAL DETAILS A. Structure models The Cu(410) surface was simulated by a five-layer p(2×1) slab (FIG. 1) with a \sim15 Å vacuum layer, where the two bottom atom layers were fixed to their bulk structure during geometry optimization. Lattice parameters of the supercell are a=7.27 Å, b=7.71 Å, and c=24.26 Å. Terrace planes of Cu(410) have a (100) surface structure and the step has a (110) facet.  FIG. 1 Surface structure of Cu(410). Miller indices of step facet and terrace plane are marked. Copper atoms at step edge are highlighted in deeper color. B. DFT calculations All DFT calculations were performed with the Vienna ab initio Simulation Package (VASP) [18]. Total energies and residual forces on each atom were converged to 10^{-4} eV and 0.02 eV/Å, respectively. A 5×4×1 k-point grid was adopted for energy and frequency calculations. For geometry optimization and transition state search, a 2×2×1 k-point grid was adopted. Transition states were located using the nudged elastic band method (NEB) [19]. All atoms except those in the two bottom layers were used to calculate the vibrational frequencies, which were used to assign experimental peaks, to confirm the transition states, and to calculate the zero-point energy (ZPE). Considering the importance of van de Waals (vdW) interaction in surface adsorption and reaction [20], the optB86b version of vdW-DF [21], which gives the best agreement with the measured adsorption energy of ethylene on Cu(100) [22], was adopted in this study. The adsorption energy (E_{\rm ads}) is defined as:  \begin{eqnarray*} E_{{\rm ads}} =E_{{\rm mol}} +E_{{\rm surf}} -E_{{\rm tot}} \end{eqnarray*} where E_{\rm tot} is the total energy of adsorbed system, E_{\rm mol} is the energy of the isolated molecule or fragments, and E_{\rm surf} means the energy of the Cu(410) surface. C. Micro-kinetic Modeling The reaction rate of each elementary step depends on the kinetic rate constant and also the surface coverage or pressure of its reactants and products:  {r_i} = {{\rm }}\left( {{k_{{{\rm f}},\mathit{i}}}\prod\limits_{j \in {{\rm IS}}} {{P_j}\prod\limits_{j \in {{\rm IS}}} {{\theta _j}} } } \right) - {{\rm }}\left( {{k_{{{\rm r}},\mathit{i}}}\prod\limits_{j \in {{\rm FS}}} {{P_j}} \prod\limits_{j \in {{\rm FS}}} {{\theta _j}} } \right) where k_{{\rm f}, i} and k_{{\rm r}, i} are the forward and backward rate constants for elementary step i respectively, P_j and \theta_j represent the partial pressure and coverage of relevant species j. Rate constants are predicted via transition state theory (TST). For example, the forward rate constant for step i is  {k_{{{\rm f}},\mathit{i}}} = A\exp \left( { - \frac{{{E_{{{\rm af}},\mathit{i}}}}}{{{k_{{\rm B}}}T}}} \right) where A={k_{\rm B}}T/h is the pre-exponential factor determined by Boltzmann's constant k_{\rm B}, Planck's constant h, and reaction temperature T, and E_{{\rm af}, i} is the calculated forward activation energy barrier for step i. In this study, ZPE correction was applied in activation energy calculations. With the reaction rate of each elementary step known, consuming rate of each species can be obtained by solving the following equation  \begin{eqnarray*} \frac{{\rm d}\theta_j}{{\rm d}t} = \sum\limits_{j \in {{\rm FS}}} {{s_n}{r_n}} - \sum\limits_{j \in {{\rm IS}}} {{s_n}{r_n}} \end{eqnarray*} where t is time and s_n is the stoichiometric factor of species j in the elementary step n. The first and the second terms in the right-hand side are from elementary steps where species j acts as a reactant and a product, respectively. Ⅲ. RESULTS A. Adsorption and the first dehydrogenation step Different adsorption configurations of ethylene on Cu(410) are investigated, which can be classified into two types: \pi-bonded (two carbon atoms bind with the same Cu atom) and di-\sigma-bonded (two carbon atoms bind with two separate Cu atoms). \pi-bonded configurations exist both at step (\pi-s) and on terrace (\pi-t). For \pi-s adsorption, there are two configurations with C-C bond parallel (\pi-s-p) and vertical (\pi-s-v) to the step. Since the terrace on Cu(410) is narrow, we also consider two configurations for the \pi-t adsorption, i.e. \pi-t-p and \pi-t-v with ethylene parallel and vertical to the step, respectively. In di-\sigma-bonded adsorption, two carbon atoms can bind with two nearest-neighboring Cu atoms (\sigma-sb) or two next-nearest-neighboring Cu atoms (\sigma-lb). For \sigma-sb, we identify three configurations with two terrace atoms (\sigma-sb-tt) or one terrace atom and one step atom (\sigma-sb-st, \sigma-sb-st') involved. As for \sigma-lb, only the configuration with two step Cu atoms involved (\sigma-lb-ss) is obtained. Optimized structures of these adsorption configurations are shown in FIG. 2 and the adsorption energies are listed in Table Ⅰ. The \pi-s-v configuration is the most stable one with an adsorption energy of 0.92 eV. The optimized C-Cu distance is 2.15 Å. di-\sigma-bonded configurations are generally less stable than \pi-bonded configurations. Our test calculations with PBE functional indicate that it significantly underestimated the adsorption energies [22].  FIG. 2 Structures of different ethylene adsorption configurations on Cu(410). (a) \pi-s-v, (b) \pi-s-p, (c) \pi-t-v, (d) \pi-t-p, (e) \sigma-sb-st, (f) \sigma-sb-tt, (g) \sigma-sb-st', (h) \sigma-lb-ss. The white, grey, and orange spheres represent hydrogen, carbon, and copper, respectively. Table Ⅰ C-Cu bond length ({d_{\rm Cu-C}}) and adsorption energy (E_{\rm ads}) for ethylene adsorption on Cu(410). Reaction energy (E_{\rm rxn}), activation energy (E_{\rm a}), and the corresponding rate constants (k) at 200 K of the first dehydrogenation step. Via the first C-H bond breaking, adsorbed ethylene can be converted to vinyl (C{\rm H}_{\rm 2}CH). Dehydrogenation energy barrier for the \pi-s-v ethylene is 1.13 eV (FIG. S1 in supplementary materials). At the transition state, the C-H distance is elongated to 1.90 Å and the shortest bond length of C-Cu is 1.97 Å. We also check the possibility of indirect dehydrogenation, i.e. the \pi-s-v configuration transfers to another configuration and then dehydrogenates. The corresponding minimum energy paths (MEP) are shown in FIG. 3(a). From \pi-s-v to \pi-s-p, the activation barrier is very low (0.03 eV), which then leads to a higher dehydrogenation barrier (1.35 eV) compared to other pathways.  FIG. 3 Minimum energy paths for indirect dehydrogenation processes via different metastable ethylene adsorption configurations. (b) Equilibrium product coverage for different indirect dehydrogenation pathways. To quantitatively evaluate these reaction pathways, an effective activation barrier is given to each pathway by fitting their time-dependent product concentration curves to that of direct dehydrogenation (FIG. S4 in supplementary materials). The results are listed in Table Ⅰ, where the pathway via \sigma-lb-ss has the lowest effective barrier (1.22 eV). At 200 K, the explicit kinetic network with all indirect dehydrogenation pathways included is constructed. The equilibrium coverage of {\rm C}_{\rm 2}$${\rm H}_{\rm 3}$ with different adsorption configurations is shown in FIG. 3(b), which clearly shows that the products generated via the $\sigma$-lb-ss pathway is several orders of magnitude larger than others. Therefore, in the following part, we only consider the direct dehydrogenation from $\pi$-s-v (${\rm C}_{\rm 2}$${\rm H}_{\rm 4}) and the indirect dehydrogenation via \sigma-lb-ss ({\rm C}_{\rm 2}$${\rm H}_{\rm 4}$$') in the first dehydrogenation step. These two processes actually have the same vinyl product structure. B. Further dehydrogenation steps and the kinetic network Further dehydrogenation involves more species on the surface (FIG. 4), including vinyl (C{\rm H}_{\rm 2}CH), vinyldene (CC{\rm H}_{\rm 2}), acetylene (CHCH), acetylidene (CCH), and dicarbon (CC). Vinyl from the first dehydrogenation step has an adsorption energy of 2.89 eV. From vinyl, both acetylene and vinylidene can be the possible dehydrogenation products. Acetylene adsorbs on a di-\sigma long-bridge step site with an adsorption energy of 1.47 eV. Vinylidene adsorbs at the step edge with an adsorption energy of 3.49 eV. Three possible adsorption configurations are found for the next dehydrogenation product, acetylidene. Two of them adopt a di-\sigma-bonded configuration, with an adsorption energy of 4.70 eV for the one parallel to the step and 4.74 eV vertical to the step. The \pi-bonded configuration at a step site has an adsorption energy of 4.86 eV. The final dehydrogenation product, dicarbon, can adopt either a long-bridge configuration with an adsorption energy of 7.28 eV or a short-bridge configuration with an adsorption energy of 6.87 eV.  FIG. 4 Possible elementary reaction steps for ethylene dehydrogenation on Cu(410) surface. White, grey, and orange spheres represent hydrogen, carbon, and copper, respectively. With these relevant species, possible elementary steps for ethylene dehydrogenation are marked in FIG. 4 and summarized in Table Ⅱ. Reaction energies and barriers are corrected by ZPE calculations. To understand the dehydrogenation mechanism, we perform a micro-kinetic analysis based on these elementary steps. Reaction rate constants are obtained by simply using a typical pre-exponential factor (4.2×10^{12} {\rm s}^{-1}). The reaction temperature is set to an experimental relevant value of 200 K [14, 15]. Since there was no hydrogen flowed into the reaction chamber in experiment [14, 15] and the hydrogen adsorption energy is low, we assume that hydrogen can desorb from the surface immediately as soon as it is generated and it will not adsorb again. Table Ⅱ Activation energy for both forward (E_{\rm af}) and reverse (E_{\rm ar}) reactions. The shortest C-Cu distance ({d_{\rm C-Cu}}) and the relevant C-H distance ({d_{\rm C-H}}) at the transition state are also listed. By supposing the initial coverage is 0.01 for ethylene and zero for all other species, we can solve the rate equations and obtain the coverage of all species as a function of time. As shown in FIG. 5, besides the reactant {\rm C}_{\rm 2}$${\rm H}_{\rm 4}$ and the product CC, the coverage of C${\rm H}_{\rm 2}$CH is the highest among all intermediate species. At the same time, from coverage change of relevant species we can obtain the reaction rate of each elementary step. Then, by integrating reaction rates, we can know how much each elementary step has happened from the beginning to the equilibrium. The results are shown in FIG. 6(a), which indicates that ${\rm C}_{\rm 2}$${\rm H}_{\rm 4}$$\rightarrow$C${\rm H}_{\rm 2}$CH$\rightarrow$CHCH$\rightarrow$CCH$'$$\rightarrowCCH''$$\rightarrow$CCH$\rightarrow$CC is the dominant pathway for ethylene dehydrogenation. The energy profile of this pathway is shown in FIG. 6(b). This conclusion is different from the previous assumption that ethylene dehydrogenate on Cu(410) should be through the long-bridge state ($\sigma$-lb-ss) [16, 17]. Notice that the unimportance of the long-bridge intermediate state is consistent with the extremely low barrier from this state to the more stable $\pi$-s-v configuration.

 FIG. 5 Coverages of all surface species as a function of time during ethylene dehydrogenation on Cu(410).
 FIG. 6 (a) Relative importance of each elementary step in the kinetic network. (b) Energy profile for the dominant kinetic pathway of ethylene dehydrogenate on Cu(410).
C. Vibrational analysis

HREELS and IRAS have been used in experiment to characterize surface species. Different vibrational peaks were obtained at different temperatures and surface coverages. For IRAS at 0.5 L@93 K [17] and HREELS at 1 L@145 K [15], experimental peaks can be largely understood with $\pi$-bonded ethylene adsorption. Two new peaks at 1330 and 2903 c${\rm m}^{-1}$ in HREELS at 51 L@193 K [15] were previously assigned to di-$\sigma$-bonded ethylene. For HREELS at 116 L@300 K [15], only a peak at 388 c${\rm m}^{-1}$ was observed and this peak persisted up to 500 K.

Since our kinetic analysis indicate that they are the most probable surface species, vibrational frequencies of adsorbed ethylene, C${\rm H}_{\rm 2}$CH and C$-$C are calculated (Table S1 in supplementary materials) to compare with experimental observations. As listed in Table Ⅲ, the vibrational frequencies of $\pi$-bonded ethylene agree well with the peaks in IRAS at 0.5 L@93 K and HREELS at 1.0 L/145 K, which confirms that ethylene adopts the $\pi$-bonded configurations at step sites at low coverage and low temperature. As for HREELS at 51 L/193 K, di-$\sigma$-bonded ethylene adsorption seems to be an explanation of the new experimental peaks. However, it is difficult to understand why the less stable di-$\sigma$-bonded configuration should be observed at higher temperature. At the same time, vibrational frequencies of C${\rm H}_{\rm 2}$CH cannot match the experimental frequencies very well in this case (1212 vs. 1123 c${\rm m}^{-1}$ and 1253 vs. 1330 c${\rm m}^{-1}$). Another possibility is ethylene and C$-$C co-adsorption (FIG. S7 in supplementary materials), where a four-member ring is formed with an ethylene adsorption energy of 0.85 eV. The calculated vibrational frequencies of this configuration at 2934, 1396, 1120, 950, and 905 c${\rm m}^{-1}$ agree well with the HREELS observation at 193 K with 51 L exposure.

Table Ⅲ Calculated vibrational modes of five adsorption states of ethylene and dicarbon compared with HREELS and IRAS measurements: 145 K HREELS 1 L, 193 K HREELS 51 L, and 93 K IRAS 0.5 L.

The vibrational frequency of C$-$C is calculated to be 346 c${\rm m}^{-1}$, which corresponds to the only peak (388 c${\rm m}^{-1}$) in HREELS at 300 K and indicates that ${\rm C}_{\rm 2}$${\rm H}_{\rm 4}$ can be dehydrogenated on Cu(410) at this temperature. At the same time, since ethylene and C$-$C co-adsorption is used to explain HREELS at 51 L/193 K, our vibrational analysis suggests that ethylene dehydrogenation can also proceed at temperature as low as 193 K.

Ⅳ. DISCUSSION AND CONCLUSIONS

Dehydrogenation of methane on Cu(410) is also investigated. The energy barriers of the four successive dehydrogenation steps are calculated as 1.20, 1.27, 0.75, and 1.41 eV, respectively. The rate-limiting step is the dehydrogenation of CH, with a barrier higher than that of ethylene decomposition (1.13 eV). The overall reaction barrier for methane dehydrogenation (2.29 eV) is also much higher than that of ethylene. The methane decomposition process is highly endothermic (1.39 eV), while the energy increase upon ethylene decomposition is moderate (0.20 or 0.31 eV). Therefore, higher temperature will be required for methane dehydrogenation on Cu(410).

In summary, on the basis of vdW corrected DFT calculations, the adsorption and dehydrogenation of ethylene on Cu(410) have been systematically investigated. The most stable configuration of ethylene adsorption is $\pi$-s-v, which is also confirmed by vibrational analysis. Micro-kinetic analysis suggests that ethylene prefers to dehydrogenate directly from the most stable configuration rather than proceed through the long-bridge metastable state. Ethylene adsorption on C$-$C may contribute to the peaks in HREELS at 193 K, which implies that dehydrogenation proceeds at such low temperature.

Supplementary materials: More detailed information on optimized geometries, microkinetic models, and results for C${\rm H}_{\rm 4}$ is given.

Ⅴ. ACKNOWLEDGMENTS

This work is partially supported by the National Natural Science Foundation of China (No.21473167 and No.21173202) and the National Key Research and Development Program of China (No.2016YFA0200600), the Fundamental Research Funds for the Central Universities (WK3430000005), and China Scholarship Council (No.201706345015). Computational resources of Super-computing Center of University of Science and Technology of China, Guangzhou and Shanghai Supercomputer Centers are also acknowledged.

Reference
 [1] A. S. Crampton, M. D. Rötzer, F. F. Schweinberger, B. Yoon, U. Landman, and U. Heiz, J. Catal. 333 , 51 (2016). DOI:10.1016/j.jcat.2015.10.023 [2] G. T. K. K. Gunasooriya, E. G. Seebauer, and M. Saeys, ACS Catal. 7 , 1966 (2017). DOI:10.1021/acscatal.6b02906 [3] K. Shimamura, Y. Shibuta, S. Ohmura, R. Arifin, and F. Shimojo, J. Phys.:Condens. Matter. 28 , 145001 (2016). DOI:10.1088/0953-8984/28/14/145001 [4] J. Andersin, N. Lopez, and K. Honkala, J. Phys. Chem. C 113 , 8278 (2009). [5] K. Kousi, N. Chourdakis, H. Matralis, D. Kontarides, C. Papadopoulou, and X. Verykios, Appl. Catal. A 518 , 129 (2016). DOI:10.1016/j.apcata.2015.11.047 [6] H. S. Bengaard, J. K. Nørskov, J. Sehested, B. S. Clausen, L. P. Nielsen, A. M. Molenbroek, and J. R. Rostrup-Nielsen, J. Catal. 209 , 365 (2002). DOI:10.1006/jcat.2002.3579 [7] H. Tetlow, J. P. de Boer, I. J. Ford, D. D. Vvedensky, J. Coraux, and L. Kantorovich, Phys. Rep. 542 , 195 (2014). DOI:10.1016/j.physrep.2014.03.003 [8] P. Wu, W. H. Zhang, Z. Y. Li, and J. L. Yang, Small 10 , 2114 (2014). DOI:10.1002/smll.201470064 [9] Z. Y. Qiu, P. Li, Z. Y. Li, and J. L. Yang, Acc. Chem. Res. 51 , 728 (2018). DOI:10.1021/acs.accounts.7b00592 [10] D. Arvanitis, L. Wenzel, and K. Baberschke, Phys. Rev. Lett. 59 , 2435 (1987). DOI:10.1103/PhysRevLett.59.2435 [11] H. Öström, A. Föhlisch, M. Nyberg, M. Weinelt, C. Heske, L. G. M. Pettersson, and A. Nilsson, Surf. Sci. 559 , 85 (2004). DOI:10.1016/j.susc.2004.04.041 [12] R. Linke, C. Becker, T. Pelster, M. Tanemura, and K. Wandelt, Surf. Sci. 377 , 655 (1997). [13] D. Yamazaki, M. Okada, F. C. Franco Jr, and T. Kasai, Surf. Sci. 605 , 934 (2011). DOI:10.1016/j.susc.2011.02.010 [14] T. Kravchuk, L. Vattuone, L. Burkholder, W. T. Tysoe, and M. Rocca, J. Am. Chem. Soc. 130 , 12552 (2008). DOI:10.1021/ja802105z [15] T. Kravchuk, V. Venugopal, L. Vattuone, L. Burkholder, W. T. Tysoe, M. Smerieri, and M. Rocca, J. Phys. Chem. C 113 , 20881 (2009). DOI:10.1021/jp904794n [16] V. Venugopal, L. Vattuone, T. Kravchuk, M. Smerieri, L. Savio, J. Jupille, and M. Rocca, J. Phys. Chem. C 113 , 20875 (2009). DOI:10.1021/jp9047924 [17] T. Makino, M. Okada, and A. Kokalj, J. Phys. Chem. C 118 , 27436 (2014). DOI:10.1021/jp509228v [18] G. Kresse, and J. Furthmller, Phys. Rev. B 54 , 11169 (1996). DOI:10.1103/PhysRevB.54.11169 [19] H. Jónsson, G. Mills, and K. W. Jacobsen, Classical and Quantum Dynamics in Condensed Phase Simulations, B. J. Berne, G. Ciccotti, and D. F. Coker Eds., Singapore River Edge, NJ: World Scientific, 385(1998). [20] J. H. Choi, Z. C. Li, P. Cui, X. D. Fan, H. Zhang, C. G. Zeng, and Z. Y. Zhang, Sci. Rep. 3 , 1925 (2013). DOI:10.1038/srep01925 [21] J. Klimeš, D. R. Bowler, and A. Michaelides, Phys. Rev. B 83 , 195131 (2011). DOI:10.1103/PhysRevB.83.195131 [22] F. Hanke, M. S. Dyer, J. Björk, and M. Persson, J. Phys.:Condens. Matter. 24 , 424217 (2012). DOI:10.1088/0953-8984/24/42/424217

a. 中国科学技术大学合肥微尺度物质科学国家研究中心, 合肥 230026;
b. 中国科学技术大学中科院能量转换材料重点实验室, 合肥 230026;
c. 中国科学技术大学量子信息与量子科技前沿协同创新中心, 合肥 230026;
d. 澳洲国立大学物理与工程研究学院应用数学系, 澳大利亚堪培拉 2600