# Fileset

[1.5018065.pdf](https://mdr.nims.go.jp/filesets/0f9e2b42-825c-4dd0-9cd2-82d0a909a2ce/download)

## Creator

[Hou, Zhufeng](https://orcid.org/0000-0002-0069-5573), [Dieb, Thaer M.](https://orcid.org/0000-0002-8111-2009), [Tsuda, Koji](https://orcid.org/0000-0002-4288-1606)

## Rights



## Other metadata

[Structure prediction of boron-doped graphene by machine learning](https://mdr.nims.go.jp/datasets/4f5b2c66-7ded-44c6-bea0-638fd212e6c2)

## Fulltext

Structure prediction of boron-doped graphene by machine learningJ. Chem. Phys. 148, 241716 (2018); https://doi.org/10.1063/1.5018065 148, 241716© 2018 Author(s).Structure prediction of boron-dopedgraphene by machine learningCite as: J. Chem. Phys. 148, 241716 (2018); https://doi.org/10.1063/1.5018065Submitted: 02 December 2017 . Accepted: 28 February 2018 . Published Online: 19 March 2018Thaer M. Dieb, Zhufeng Hou , and Koji Tsuda COLLECTIONSPaper published as part of the special topic on Data-Enabled Theoretical ChemistryARTICLES YOU MAY BE INTERESTED INSchNet – A deep learning architecture for molecules and materialsThe Journal of Chemical Physics 148, 241722 (2018); https://doi.org/10.1063/1.5019779Constructing first-principles phase diagrams of amorphous LixSi using machine-learning-assisted sampling with an evolutionary algorithmThe Journal of Chemical Physics 148, 241711 (2018); https://doi.org/10.1063/1.5017661Alchemical and structural distribution based representation for universal quantum machinelearningThe Journal of Chemical Physics 148, 241717 (2018); https://doi.org/10.1063/1.5020710https://images.scitation.org/redirect.spark?MID=176720&plid=1085727&setID=378408&channelID=0&CID=358608&banID=519992853&PID=0&textadID=0&tc=1&type=tclick&mt=1&hc=a6e1cecbc242d3b912549e1a9893d52b6202f329&location=https://doi.org/10.1063/1.5018065https://doi.org/10.1063/1.5018065https://aip.scitation.org/author/M+Dieb%2C+Thaerhttps://aip.scitation.org/author/Hou%2C+Zhufenghttp://orcid.org/0000-0002-0069-5573https://aip.scitation.org/author/Tsuda%2C+Kojihttp://orcid.org/0000-0002-4288-1606/topic/special-collections/detc2018?SeriesKey=jcphttps://doi.org/10.1063/1.5018065https://aip.scitation.org/action/showCitFormats?type=show&doi=10.1063/1.5018065http://crossmark.crossref.org/dialog/?doi=10.1063%2F1.5018065&domain=aip.scitation.org&date_stamp=2018-03-19https://aip.scitation.org/doi/10.1063/1.5019779https://doi.org/10.1063/1.5019779https://aip.scitation.org/doi/10.1063/1.5017661https://aip.scitation.org/doi/10.1063/1.5017661https://doi.org/10.1063/1.5017661https://aip.scitation.org/doi/10.1063/1.5020710https://aip.scitation.org/doi/10.1063/1.5020710https://doi.org/10.1063/1.5020710THE JOURNAL OF CHEMICAL PHYSICS 148, 241716 (2018)Structure prediction of boron-doped graphene by machine learningThaer M. Dieb,1,2,3,a) Zhufeng Hou,2,a) and Koji Tsuda1,2,3,b)1Graduate School of Frontier Sciences, The University of Tokyo, Kashiwa, Japan2National Institute for Materials Science (NIMS), Tsukuba, Japan3RIKEN, Center for Advanced Intelligence Project, Tokyo, Japan(Received 2 December 2017; accepted 28 February 2018; published online 19 March 2018)Heteroatom doping has endowed graphene with manifold aspects of material properties and boostedits applications. The atomic structure determination of doped graphene is vital to understand itsmaterial properties. Motivated by the recently synthesized boron-doped graphene with relatively highconcentration, here we employ machine learning methods to search the most stable structures of dopedboron atoms in graphene, in conjunction with the atomistic simulations. From the determined stablestructures, we find that in the free-standing pristine graphene, the doped boron atoms energeticallyprefer to substitute for the carbon atoms at different sublattice sites and that the para configuration ofboron-boron pair is dominant in the cases of high boron concentrations. The boron doping can increasethe work function of graphene by 0.7 eV for a boron content higher than 3.1%. © 2018 Author(s).All article content, except where otherwise noted, is licensed under a Creative Commons Attribution(CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5018065I. INTRODUCTIONChemical composition modification of materials by incor-porating additional elements is one of the commonly usedapproaches to tune the material properties. To achieve thedesired material properties, not only the optimal composi-tion but also the spatial distribution of incorporated atomsis demanded to be discovered. Because of the complexityintroduced by the interplay between structural and chemi-cal degrees of freedom, the atomic structure search of dopantatoms in the host material is very challenging, particularly forthe non-diluted doping level.Graphene, a single layer of sp2-bonded carbon atomsarranged in a honeycomb lattice, is the most explored repre-sentative of two-dimensional (2D) materials.1,2 However, theelectrical conduction of pristine monolayer graphene cannotbe switched off because of its zero band gap in the electronicdensity of states, limiting the range of potential applications.2,3The perfectly flat graphene is chemically inert and not of prac-tical interest in the applications involving chemical reactions.4To expand the applications of graphene, many efforts havebeen made to modify graphene by chemical functionalization,chemical doping, structure engineering, and so on.2,4–6 Amongthese approaches, incorporation of heteroatoms into grapheneis shown to be a versatile method for controllable tuning of itsphysical and chemical properties.5–10 Boron (B) and nitrogen(N) have attracted much great attention for chemical dopingof graphene since they are neighboring to carbon (C) and theiratomic radii are similar to that of carbon.11,12 Compared withnitrogen, substitutional doping of boron in graphene was pre-dicted to be more energetically favorable13,14 and found toexhibit a p-type doping effect.14,15 Recent experiments havea)T. M. Dieb and Z. Hou contributed equally to this work.b)tsuda@k.u-tokyo.ac.jpshown that boron can be embedded in graphene at higher con-centrations (up to 20%).10,14,16–18 It was also found that thehexagonal structure exists in boron carbides (BxC1�x), thinfilms with B contents less than 0.5.19 For BxC1�x, thin filmswith B content at x = 0.25, an ordered layer structure was pro-posed.20,21 These make boron-doped graphene (B-graphene)a unique and very attractive material from both fundamentaland practical viewpoints. However, the local structural formof B-graphene with relatively high B concentration is notyet well understood. The spatial distribution of doped B ingraphene is strongly affected by the synthesis methods andconditions. For the epitaxial B-graphene grown by chemicalvapor deposition, the B dopants in graphene can be a pref-erential substitution of carbon in only one of the graphenesublattices10,18 and be completely random,22 depending onthe metal substrate used. To understand the energetic stabil-ity and to find out possible ordered structures of B-graphene,we carry out a global search for the local structure configu-rations of the B dopants in graphene with a series of concen-trations by means of the supercell model and the atomisticsimulations.In the case of substitutional doping of B for C in defect-free graphene, the number of possible structure configurationscan be written as Cmn+m if the symmetry is broken, where nand m are the numbers of C and B atoms, respectively. Itis evident that the exhaustive search of all possible structureconfigurations for multiple B dopants in graphene would beprohibitive from the efficiency. To search the stable atomicstructures of B-graphene more efficiently, in this work we pro-pose to use a machine learning approach that utilizes MonteCarlo tree search,23 which was recently utilized for materialsand chemical design24,25 with a rollout schema depending onBayesian optimization (BO).26 Based on the structure searchof B-graphene, we discuss the performance of this proposedscheme of global optimization.0021-9606/2018/148(24)/241716/7 148, 241716-1 © Author(s) 2018https://doi.org/10.1063/1.5018065http://creativecommons.org/licenses/by/4.0/https://doi.org/10.1063/1.5018065mailto:tsuda@k.u-tokyo.ac.jphttp://crossmark.crossref.org/dialog/?doi=10.1063/1.5018065&domain=pdf&date_stamp=2018-03-19241716-2 M. Dieb, Hou, and Tsuda J. Chem. Phys. 148, 241716 (2018)The remainder of this paper is organized as follows. InSec. II, we briefly introduce the machine learning methodsused for the search of atomic structures of B-graphene. Thedetails of computation setup are given in Sec. III. The resultsfor the optimization performance of the machine learningmethods, the stabilities of B-graphene, and the electronic struc-tures of B-graphene are presented in Sec. IV. The conclusionsare given in Sec. V.II. MACHINE LEARNING BASEDOPTIMIZATION METHODDetermination of optimal material structure with cer-tain quality metrics traditionally depended on the experi-ence of domain knowledge experts and trial-and-error exper-iments. Several machine learning-based methods have beenproposed to accelerate this process with as few experimentsas possible.24,27–31 In such methods, the structure predictionis often formulated as a selection of optimal solution froma candidates space that maximizes or minimizes a black-box function (usually the target property).32,33 Experimentsare commonly replaced by simulators such as first-principlescalculations.Currently available efficient methods such as Bayesianoptimization (BO)26,34 are not scalable enough. Due to itsexceptional performance in computer Go game,23,35 MonteCarlo tree search (MCTS)23 has recently gained attention inmaterial design24 offering a superior scalability with efficiencytrade-off.MCTS employs a shallow search tree where a node is apossible assignment of an atom into a position in the struc-ture. Iteratively, the tree expands towards the optimal solutionin 4 steps: selection, expansion, simulation, and backprop-agation. A full solution is obtained from this shallow treeusing the random rollout (completion) technique in the sim-ulation step. To increase the MCTS efficiency, we proposeto use Bayesian learning to engineer the rollout24 instead ofthe random selection. We briefly discuss this rollout oper-ation here. See the supplementary material for details onthe proposed method. The source code is freely available athttps://github.com/tsudalab/MDTS.Bayesian optimization methods maintain a surrogatemodel of the black-box function, most commonly Gaussianprocess (GP). A pool of candidate S′ is generated at a nodeof level p` where each data point represents a full structureof N positions, with positions p1, . . ., p` being determinedalong the path from the root to the selected node and posi-tions p`+1, . . ., pN being randomly generated. Data points inS′ are then vectors of binary values (0, 1). If no data pointfrom S′ is previously observed, GP starts with an initial setof randomly selected data points from S′. GP is updatedas more data points are observed. An acquisition functionis then used to determine the next optimal solution. Withincomputational budget, the selected solutions are returned(Fig. 1).III. COMPUTATIONAL DETAILSTo study B-graphene, a supercell constructed by the4 × 4 extension of the hexagonal unit cell of graphene wasemployed and the substitution of carbon by boron was con-sidered. To avoid the spurious interaction between graphenelayers, a vacuum thickness in the supercell was set to 20.0 Å.The in-plane lattice constant of graphene supercell wasfixed at 4a0, where a0 = 2.464 Å for the calculated lat-tice constant of graphene. Figure 2 shows the atomic struc-ture of the graphene supercell used in the present study.The number of B dopants was considered from 1 to 10,which correspond to the boron concentration variation from3.125% to 31.25%. The number of possible structure con-figurations for different numbers of doped B atoms in a4 × 4 graphene supercell is summarized in Table I. Becauseof the large number of structure configurations for multiple Bdopants (nB > 4) in graphene, we employed a combinatorialFIG. 1. Monte Carlo tree search(MCTS) with Bayesian rollout forbinary atom assignment. MCTS uses ashallow tree search. At a certain itera-tion, only partial solution is determinedin the tree. To obtain full solutions,a pool of full solution candidates isenumerated. GP is fed with previousobservations if available or initialrandom selection (red triangles). Anacquisition function is applied todetermine the next optimal observationfrom the pool.ftp://ftp.aip.org/epaps/journ_chem_phys/E-JCPSA6-148-017891https://github.com/tsudalab/MDTS241716-3 M. Dieb, Hou, and Tsuda J. Chem. Phys. 148, 241716 (2018)FIG. 2. The atomic structures of a 4 × 4 supercell of monolayer graphene (a)without and (b) with one doped B atom.structure-generation approach for the local-level modeling ofatomic substitutions and partial occupancies in crystals toreduce the number of potential configurations.36 Our proposedalgorithm is then applied to search the stable structure config-urations of B-doped graphene. During the structure search,the atomic positions in each structure configuration were opti-mized by the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method.37 The energy and atomic forcewere calculated by the classical force field method,38 as imple-mented in the ATK-classical simulator of Atomistix ToolKit(ATK).39,40 The interactions between the atoms in the systemunder study were described by the Tersoff bond order poten-tials developed by Matsunaga et al.41 for C–B and C–C inter-actions, which have been widely used to study the B-dopedcarbon nanostructures. For the top ten structures with lowerenergies found by the global search methods for each caseof multiple B dopants, we carried out density functional the-ory (DFT) calculations to optimize the atomic position andcalculate the energy.Density functional theory calculations have been per-formed with the Vienna Ab initio Simulation Package(VASP).42,43 Perdew-Burke-Ernzerh of (PBE) exchange-correlation functional44 within the generalized gradientapproximations (GGA) was used in conjunction with theprojector-augmented wave (PAW) method.45,46 A plane wavebasis set with an energy cutoff of 500 eV was used and ak-point mesh was sampled with a 11 × 11 × 1 Monkhorst-Pack scheme.47 The criterions for total energy and forceconvergence were set to 10�5 eV and 10�2 eV/Å, respec-tively. Because of the fairly delocalized nature for the defectstates induced by doped boron atoms,48,49 the spin polarizationof B-graphene is energetically unfavorable, which has beenconfirmed in the test calculations of spin-polarized DFT.Therefore, we report here the results from the non-spin-polarized calculations. Since the well-known underestima-tion of band gap in the GGA-PBE calculations due to self-interaction errors, we have also performed the HSE06 hybridfunctional calculations50 to study the electronic structures ofB-graphene in some cases.IV. RESULTS AND DISCUSSIONA. Obtaining optimal structuresAs listed in Table I, the number of symmetry non-equivalent atomic structures for multiple doped B atomsconsidered in a 4 × 4 graphene supercell can reach severalhundred thousands. For such large space of candidate struc-tures of B-graphene, the atomic structure search was carriedout by the MCTS with Bayesian rollout. We present in Fig. 3the evolution of lattice energy during the structure search. Inthe cases considered here, we can see that the change of latticeenergies for found optimal structures is less than 0.02 eV persupercell as the search reaches the predefined maximum steps,which indicates a good convergence for the structure search.For the cases of search space with thousand candidate struc-tures, the finding of optimal structure takes 40% of the searchspace. But for the case of search space with candidate struc-tures more than one hundred thousand, the optimal structureof B-graphene could be found within 2% of the search space.Therefore, our proposed optimization scheme might be verysuitable for the larger search space that is hard to be treated, forexample, the use of larger supercell size with a kept numberof dopant atoms.B. Structure stability of boron-doped grapheneTo evaluate the stabilities of doped B atoms in graphene,we calculate their formation energies according to the follow-ing definition:∆Eform =1mB[Et(BG) − Et(PG) − mBµB + mBµC], (1)where Et(BG) and Et(PG) are the total energies of graphenesupercell with and without B dopants, mB is the number ofdoped B atoms, and µB and µC are the chemical potentialsof boron and carbon, which are taken as the total energiesper atom of α-boron crystal and pristine monolayer graphene,respectively. The formation of single B substitution in the4 × 4 graphene supercell is about 1.186 eV, which is in goodagreement with the result (1.12 eV) reported in the previousstudy.51It is well known that the honeycomb lattice of grapheneconsists of two sublattices, namely, A and B sublattices. Tounderstand the occupancy preference of multiple-doped Batoms in graphene, it is necessary to first examine the inter-action of two B substitutions in graphene. For this purpose,we calculate the interaction energy of two B substitutions asdefined in the following equation:∆Eint = Et(2B) + Et(PG) − 2Et(1B), (2)TABLE I. The number of possible structure configurations for mB (mB = 1, . . ., 10) doped B atoms in a4 × 4 graphene supercell after merging symmetry.36 The corresponding concentration of doped B in graphene iscB = mB/32.mB 1 2 3 4 5 6 7 8 9 10cB (%) 3.125 6.25 9.375 12.5 15.625 18.75 21.875 25 28.125 31.25Nconf. 1 8 37 241 1129 5002 17 929 55 817 147 362 338 741241716-4 M. Dieb, Hou, and Tsuda J. Chem. Phys. 148, 241716 (2018)FIG. 3. The evolution of lattice energy during the structure search of B-graphene by the optimization scheme based on MCTS with Bayesian rollout for thedifferent numbers of doped B atoms in a 4 × 4 supercell: (a) mB = 5, (b) mB = 6, (c) mB = 7, (d) mB = 8, (e) mB = 9, and (f) mB = 10.where Et(PG), Et(1B), and Et(2B) are the total energies ofgraphene supercells with zero, one, and two doped B atoms,respectively. The positive (negative) sign of the interactionenergy defined in Eq. (2) indicates two substitutional B dopantsin graphene that repel (attract) each other. It also implies theincrease (decreasing) in the formation energy of two B substi-tutions with respect to that of single B substitution. In the4 × 4 graphene supercell considered here, the most stableconfiguration of two B substitutions is depicted in Fig. 4(a),where two doped B atoms occupy the different sublattices,and the calculated interaction energy is about 0.022 eV. Thisindicates that the interaction of two substitutional B atomsin graphene is repulsive. For two B substitutions at the firstnearest-neighboring (NN) lattice sites (called ortho configura-tion) in graphene, they show the strongest repulsive interactionand the corresponding interaction energy is 1.262 eV. For twoB substitutions at the second and third NN lattice sites (calledmeta and para configurations, respectively) in graphene, theirinteraction energies are 0.507 eV and 0.075 eV, respectively.As the para configuration of the B-B pair is more stablethan the meta configuration, it has been found in the chem-ically synthesized B-doped nanographene by using the boron-containing polycyclic aromatic hydrocarbon (PAH).52,53 TheB-graphene with such less stable meta configuration of B-Bpair was synthesized by chemical vapor deposition (CVD) ona nickel or cobalt substrate that has strong interaction withgraphene.10,18The atomic structures for the most stable configurationsof B-graphene at different B concentrations are depictedin Fig. 4. They have been confirmed in the calculationsthat take into account the relaxation of in-plane lattice con-stants. We have verified their stabilities in the calculations of241716-5 M. Dieb, Hou, and Tsuda J. Chem. Phys. 148, 241716 (2018)FIG. 4. The most stable structure con-figurations for doped B atoms in the4 × 4 supercell of graphene at differentconcentrations. From (a) to (i), the num-ber of doped B atoms changes from twoto ten. The gray and red balls representC and B atoms, respectively.local-density approximations (LDAs) and meta-GGA. Thedetailed comparison of the results is given in the supple-mentary material. The corresponding formation energies ofB substitutions are presented in Fig. 5. In these structure con-figurations, the para configuration of the B-B pair is dominant,indicating the substitutional B atoms prefer to occupy dif-ferent sublattices in the free standing pristine graphene. Theunbalanced sublattice doping observed in experiments maybe attributed to the strong interaction of metal substrate andgraphene.10,18 In particular, it is noticed that the formationenergy of four B substitutions is lower than those of three andfive ones. In the case of eight B substitutions, which corre-sponds to BC3, its formation energy is lower than those ofseven and nine ones. For the most stable structure configu-ration of eight B substitutions, as shown in Fig. 4(g), all ofthe benzene-like C6 rings are separated by B-B pairs and thusits aromaticity is very high. This structure configuration hasbeen proposed in the literature20,21 for the synthesized BC3compound.FIG. 5. The formation energy of doped B atom in a 4×4 supercell of grapheneas a function of the number of doped B atoms.C. Electronic structures of boron-doped grapheneTo explore the electronic structure of graphene modifiedby B doping, we present in Fig. 6 the total density of states(DOS) of B-graphenes at a series of B concentrations. It is wellknown that the Fermi level (EF) of perfect graphene coincideswith the Dirac point. Because boron has one less electron thancarbon, the B substitution leads to the shifting of EF down intothe valence bands. We also notice that the B doping opens upa band gap at the Dirac point of graphene, which is supportedby the electrical conductivity measurement of B-graphene.14In the case of eight B atoms doped in the 4× 4 graphene super-cell, a band gap is open around the EF and thus this structureof B-graphene is a semiconductor, which is consistent withthe previous results.54,55 The band gaps of BC3 predicted bythe GGA-PBE and HSE06 hybrid functional calculations inthe present study are about 0.17 eV and 1.38 eV, respectively.To quantitatively characterize the doping level of carriers, wepresent the work function of B-graphene in Fig. 7. The workfunction of perfect monolayer graphene predicted by the GGA-PBE calculations in the present study is about 4.25 eV, whichis in good agreement with the experiment results (4.272 eV forCVD-grown graphene measured by UV photoelectron spec-troscopy56 and 4.57 ± 0.05 eV for mechanically exfoliatedgraphene measured by a scanning Kelvin probe microscope,57respectively.). It also agrees well with the previous result(4.20 eV) of GGA-PBE calculations in the study of Lazaret al.58 We can see that the boron doping increases the workfunction of graphene significantly. The similar trend was alsofound in the study of Lazar et al.58 Considering the underes-timation of band gap in the GGA-PBE calculations, we havealso performed the HSE06 hybrid functional calculations tocheck the work functions of perfect monolayer graphene andB-graphene, which includes the cases of singe B substitutionin 4 × 4 supercell and BC3. The obtained work functions forthese three cases are 4.23, 5.23, and 5.80 eV, respectively.ftp://ftp.aip.org/epaps/journ_chem_phys/E-JCPSA6-148-017891ftp://ftp.aip.org/epaps/journ_chem_phys/E-JCPSA6-148-017891241716-6 M. Dieb, Hou, and Tsuda J. Chem. Phys. 148, 241716 (2018)FIG. 6. The total density of states(DOS) of B-graphene with differentnumbers of doped B atoms in a 4 × 4supercell for their atomic structuresshown in Figs. 2 and 4. The Fermi levelis set as zero and indicated by the dashedline.FIG. 7. The work function of B-graphene as a function of the number ofdoped B atoms in a 4 × 4 supercell. In the case of two B dopants, the workfunction for the metastable structure of B-B pair in the para configuration isalso given.Depending on the doped B concentration, the GGA-PBE cal-culations show that the work function of graphene can be tunedby 0.7 eV–1.0 eV.V. CONCLUSIONIn summary, we have employed the Monte Carlo treesearch (MCTS) with Bayesian rollout to search the stablestructures of B-graphene for the boron concentration up to31.25%. Compared with the sole use of Bayesian optimiza-tion, this integrated optimization method shows a superiorityof better scalability. Our results show that in the free-standinggraphene, the doped boron atoms energetically prefer to sub-stitute for the carbon atoms at different sublattices. For B-graphene with high boron concentration, the para B-B pairsare dominant. Because of the repulsive interaction betweenboron substitutions in graphene, the doped boron would tendto be segregated. The doped boron can open a band gap atthe Dirac point of graphene. Particularly in the concentrationof doped B at 25%, namely, BC3, the found stable structureexhibits semiconducting behavior. We also find that boron dop-ing can lead to an increase of the work function of graphenein the range between 0.7 eV and 1.0 eV. Our results would bevery helpful to further explore the application of B-graphene.In addition, the proposed method can be applied to multipleatom types assignment problems such as boron and nitro-gen codoped graphene by allowing more than two branchesin the search tree and using a descriptor based on one-hotencoding.SUPPLEMENTARY MATERIALSee supplementary material for the details on the pro-posed machine learning method and the detailed comparisonof the results obtained by the different exchange-correlationfunctionals in DFT calculations.ACKNOWLEDGMENTSThis work was partially supported by “Materials researchby Information Integration” Initiative (MI2I) project of theSupport Program for Starting Up Innovation Hub from JapanScience and Technology Agency (JST). The computation inthis study was performed on Numerical Materials Simulatorat NIMS.The authors would like to thank Kei Terayama andShenghong Ju for their fruitful discussion. Also, we are thank-ful for anonymous referees for their comments and suggestionsto improve the manuscript.ftp://ftp.aip.org/epaps/journ_chem_phys/E-JCPSA6-148-017891241716-7 M. Dieb, Hou, and Tsuda J. Chem. Phys. 148, 241716 (2018)1A. K. Geim and K. S. Novoselov, Nat. Mater. 6, 183 (2007).2M. F. Craciun, I. Khrapach, M. D. Barnes, and S. Russo, J. Phys.: Condens.Matter 25, 423201 (2013).3P. Avouris, Nano Lett. 10, 4285 (2010).4A. Eftekhari and H. Garcia, Mater. Today Chem. 4, 1 (2017).5L. K. Putri, W.-J. Ong, W. S. Chang, and S.-P. Chai, Appl. Surf. Sci. 358, 2(2015).6H. Liu, Y. Liu, and D. Zhu, J. Mater. Chem. 21, 3335 (2011).7X. Wang, G. Sun, P. Routh, D.-H. Kim, W. Huang, and P. Chen, Chem. Soc.Rev. 43, 7067 (2014).8J. Duan, S. Chen, M. Jaroniec, and S. Z. Qiao, ACS Catal. 5, 5207 (2015).9H. Cui, Z. Zhou, and D. Jia, Mater. Horiz. 4, 7 (2017).10D. Y. Usachov, A. V. Fedorov, A. E. Petukhov, O. Y. Vilkov, A. G. Rybkin,M. M. Otrokov, A. Arnau, E. V. Chulkov, L. V. Yashina, M. Farjam,V. K. Adamchuk, B. V. Senkovskiy, C. Laubschat, and D. V. Vyalikh, ACSNano 9, 7314 (2015).11H. Wang, T. Maiyalagan, and X. Wang, ACS Catal. 2, 781 (2012).12S. Agnoli and M. Favaro, J. Mater. Chem. A 4, 5002 (2016).13L. S. Panchakarla, K. S. Subrahmanyam, S. K. Saha, A. Govindaraj,H. R. Krishnamurthy, U. V. Waghmare, and C. N. R. Rao, Adv. Mater.21, 4726 (2009).14Y.-B. Tang, L.-C. Yin, Y. Yang, X.-H. Bo, Y.-L. Cao, H.-E. Wang,W.-J. Zhang, I. Bello, S.-T. Lee, H.-M. Cheng, and C.-S. Lee, ACS Nano 6,1970 (2012).15W.-C. Yen, H. Medina, J.-S. Huang, C.-C. Lai, Y.-C. Shih, S.-M. Lin,J.-G. Li, Z. M. Wang, and Y.-L. Chueh, J. Phys. Chem. C 118, 25089(2014).16W. Zhao, J. Gebhardt, K. Gotterbarm, O. Höfert, C. Gleichweit, C. Papp,A. Görling, and H.-P. Steinrück, J. Phys.: Condens. Matter 25, 445002(2013).17J. Gebhardt, R. J. Koch, W. Zhao, O. Höfert, K. Gotterbarm, S. Mammadov,C. Papp, A. Görling, H.-P. Steinrück, and T. Seyller, Phys. Rev. B 87, 155437(2013).18D. Y. Usachov, A. V. Fedorov, O. Y. Vilkov, A. E. Petukhov, A. G. Rybkin,A. Ernst, M. M. Otrokov, E. V. Chulkov, I. I. Ogorodnikov, M. V. Kuznetsov,L. V. Yashina, E. Y. Kataev, A. V. Erofeevskaya, V. Y. Voroshnin, V. K.Adamchuk, C. Laubschat, and D. V. Vyalikh, Nano Lett. 16, 4535(2016).19I. Caretti, R. Gago, J. M. Albella, and I. Jiménez, Phys. Rev. B 77, 174109(2008).20J. Kouvetakis, R. B. Kaner, M. L. Sattler, and N. Bartlett, J. Chem. Soc.,Chem. Commun. 1986, 1758.21H. Yanagisawa, T. Tanaka, Y. Ishida, M. Matsue, E. Rokuta, S. Otani, andC. Oshima, Phys. Rev. Lett. 93, 177003 (2004).22L. Zhao, M. Levendorf, S. Goncher, T. Schiros, L. Pálová, A.Zabet-Khosousi, K. T. Rim, C. Gutiérrez, D. Nordlund, C. Jaye, M.Hybertsen, D. Reichman, G. W. Flynn, J. Park, and A. N. Pasupathy, NanoLett. 13, 4659 (2013).23C. B. Browne, E. Powley, D. Whitehouse, S. M. Lucas, P. I. Cowling,P. Rohlfshagen, S. Tavener, D. Perez, S. Samothrakis, and S. Colton, IEEETrans. Comput. Intell. AI Games 4, 1 (2012).24T. M. Dieb, S. Ju, K. Yoshizoe, Z. Hou, J. Shiomi, and K. Tsuda, Sci.Technol. Adv. Mater. 18, 498 (2017).25X. Yang, J. Zhang, K. Yoshizoe, K. Terayama, and K. Tsuda, Sci. Technol.Adv. Mater. 18, 972 (2017).26D. R. Jones, M. Schonlau, and W. J. Welch, J. Global Optim. 13, 455 (1998).27D. Reker and G. Schneider, Drug Discov. Today 20, 458 (2015).28A. R. Oganov and C. W. Glass, J. Chem. Phys. 124, 244704 (2006).29M. Ahmadi, M. Vogt, P. Iyer, J. Bajorath, and H. Fr öhlich, J. Chem. Inf.Model. 53, 553 (2013).30S. Ju, T. Shiga, L. Feng, Z. Hou, K. Tsuda, and J. Shiomi, Phys. Rev. X 7,021024 (2017).31R. Gómez-Bombarelli, D. K. Duvenaud, J. M. Hernández-Lobato,J. Aguilera-Iparraguirre, T. D. Hirzel, R. P. Adams, and A. Aspuru-Guzik,ACS Cent. Sci. 4(2), 268 (2018).32A. Seko, A. Togo, H. Hayashi, K. Tsuda, L. Chaput, and I. Tanaka, Phys.Rev. Lett. 115, 205901 (2015).33P. V. Balachandran, D. Xue, J. Theiler, J. Hogden, and T. Lookman, Sci.Rep. 6, 19660 (2016).34J. Snoek, H. Larochelle, and R. P. Adams, in Advances in Neural InformationProcessing Systems 25, edited by F. Pereira, C. J. C. Burges, L. Bottou, andK. Q. Weinberger (Curran Associates, Inc., 2012), pp. 2951–2959.35D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. van den Driessche,J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot, S.Dieleman, D. Grewe, J. Nham, N. Kalchbrenner, I. Sutskever, T. Lillicrap,M. Leach, K. Kavukcuoglu, T. Graepel, and D. Hassabis, Nature 529, 484(2016).36K. Okhotnikov, T. Charpentier, and S. Cadars, J. Cheminf. 8, 17 (2016).37J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed. (Springer,New York, 2006).38J. Schneider, J. Hamaekers, S. T. Chill, S. Smidstrup, J. Bulin, R. Thesen,A. Blom, and K. Stokbro, Modell. Simul. Mater. Sci. Eng. 25, 085007(2017).39See www.quantumwise.com for QuantumWise A/S, Atomistix ToolKitversion 2016.01.40M. Griebel and J. Hamaekers, Comput. Methods Appl. Mech. Eng. 193,1773 (2004).41K. Matsunaga, C. Fisher, and H. Matsubara, Jpn. J. Appl. Phys., Part 2 39,L48 (2000).42G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996).43G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).44J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).45P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).46G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).47H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).48K. E. Kweon and G. S. Hwang, Phys. Rev. B 82, 195439 (2010).49B. Zheng, P. Hermet, and L. Henrard, ACS Nano 4, 4165 (2010).50J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 124, 219906(2006).51I. Choudhuri, N. Patra, A. Mahata, R. Ahuja, and B. Pathak, J. Phys. Chem.C 119, 24827 (2015).52C. Dou, S. Saito, K. Matsuo, I. Hisaki, and S. Yamaguchi, Angew. Chem.,Int. Ed. 51, 12206 (2012).53S. Osumi, S. Saito, C. Dou, K. Matsuo, K. Kume, H. Yoshikawa, K. Awaga,and S. Yamaguchi, Chem. Sci. 7, 219 (2016).54D. Tomanek, R. M. Wentzcovitch, S. G. Louie, and M. L. Cohen, Phys. Rev.B 37, 3134 (1988).55X. Chen, Y. Yao, H. Yao, F. Yang, and J. Ni, Phys. Rev. B 92, 174503(2015).56S. Bae, H. Kim, Y. Lee, X. Xu, J.-S. Park, Y. Zheng, J. Balakrishnan, T. Lei,H. Ri Kim, Y. I. Song, Y.-J. Kim, K. S. Kim, B. Özyilmaz, J.-H. Ahn,B. H. Hong, and S. Iijima, Nat. Nanotechnol. 5, 574 (2010).57Y.-J. Yu, Y. Zhao, S. Ryu, L. E. Brus, K. S. Kim, and P. Kim, Nano Lett. 9,3430 (2009).58P. Lazar, R. Zboril, M. Pumera, and M. Otyepka, Phys. Chem. Chem. Phys.16, 14231 (2014).https://doi.org/10.1038/nmat1849https://doi.org/10.1088/0953-8984/25/42/423201https://doi.org/10.1088/0953-8984/25/42/423201https://doi.org/10.1021/nl102824hhttps://doi.org/10.1016/j.mtchem.2017.02.003https://doi.org/10.1016/j.apsusc.2015.08.177https://doi.org/10.1039/c0jm02922jhttps://doi.org/10.1039/c4cs00141ahttps://doi.org/10.1039/c4cs00141ahttps://doi.org/10.1021/acscatal.5b00991https://doi.org/10.1039/c6mh00358chttps://doi.org/10.1021/acsnano.5b02322https://doi.org/10.1021/acsnano.5b02322https://doi.org/10.1021/cs200652yhttps://doi.org/10.1039/c5ta10599dhttps://doi.org/10.1002/adma.200901285https://doi.org/10.1021/nn3005262https://doi.org/10.1021/jp508365hhttps://doi.org/10.1088/0953-8984/25/44/445002https://doi.org/10.1103/physrevb.87.155437https://doi.org/10.1021/acs.nanolett.6b01795https://doi.org/10.1103/physrevb.77.174109https://doi.org/10.1039/c39860001758https://doi.org/10.1039/c39860001758https://doi.org/10.1103/physrevlett.93.177003https://doi.org/10.1021/nl401781dhttps://doi.org/10.1021/nl401781dhttps://doi.org/10.1109/tciaig.2012.2186810https://doi.org/10.1109/tciaig.2012.2186810https://doi.org/10.1080/14686996.2017.1344083https://doi.org/10.1080/14686996.2017.1344083https://doi.org/10.1080/14686996.2017.1401424https://doi.org/10.1080/14686996.2017.1401424https://doi.org/10.1023/a:1008306431147https://doi.org/10.1016/j.drudis.2014.12.004https://doi.org/10.1063/1.2210932https://doi.org/10.1021/ci3004682https://doi.org/10.1021/ci3004682https://doi.org/10.1103/physrevx.7.021024https://doi.org/10.1021/acscentsci.7b00572https://doi.org/10.1103/physrevlett.115.205901https://doi.org/10.1103/physrevlett.115.205901https://doi.org/10.1038/srep19660https://doi.org/10.1038/srep19660https://doi.org/10.1038/nature16961https://doi.org/10.1186/s13321-016-0129-3https://doi.org/10.1088/1361-651x/aa8ff0http://www.quantumwise.comhttps://doi.org/10.1016/j.cma.2003.12.025https://doi.org/10.1143/jjap.39.l48https://doi.org/10.1016/0927-0256(96)00008-0https://doi.org/10.1103/physrevb.54.11169https://doi.org/10.1103/physrevlett.77.3865https://doi.org/10.1103/physrevb.50.17953https://doi.org/10.1103/physrevb.59.1758https://doi.org/10.1103/physrevb.13.5188https://doi.org/10.1103/physrevb.82.195439https://doi.org/10.1021/nn1002425https://doi.org/10.1063/1.2204597https://doi.org/10.1021/acs.jpcc.5b07359https://doi.org/10.1021/acs.jpcc.5b07359https://doi.org/10.1002/anie.201206699https://doi.org/10.1002/anie.201206699https://doi.org/10.1039/c5sc02246khttps://doi.org/10.1103/physrevb.37.3134https://doi.org/10.1103/physrevb.37.3134https://doi.org/10.1103/physrevb.92.174503https://doi.org/10.1038/nnano.2010.132https://doi.org/10.1021/nl901572ahttps://doi.org/10.1039/c4cp01638f