Structure prediction of boron-doped graphene by machine learning Open Access

Chemical composition modification of materials by incorporating additional elements is one of the commonly used approaches to tune the material properties. To achieve the desired material properties, not only the optimal composition but also the spatial distribution of incorporated atoms is demanded to be discovered. Because of the complexity introduced by the interplay between structural and chemical degrees of freedom, the atomic structure search of dopant atoms in the host material is very challenging, particularly for the non-diluted doping level.

Graphene, a single layer of sp²-bonded carbon atoms arranged in a honeycomb lattice, is the most explored representative of two-dimensional (2D) materials.^1,2 However, the electrical conduction of pristine monolayer graphene cannot be switched off because of its zero band gap in the electronic density of states, limiting the range of potential applications.^2,3 The perfectly flat graphene is chemically inert and not of practical interest in the applications involving chemical reactions.⁴ To expand the applications of graphene, many efforts have been made to modify graphene by chemical functionalization, chemical doping, structure engineering, and so on.^2,4–6 Among these approaches, incorporation of heteroatoms into graphene is shown to be a versatile method for controllable tuning of its physical and chemical properties.^5–10 Boron (B) and nitrogen (N) have attracted much great attention for chemical doping of graphene since they are neighboring to carbon (C) and their atomic radii are similar to that of carbon.^11,12 Compared with nitrogen, substitutional doping of boron in graphene was predicted to be more energetically favorable^13,14 and found to exhibit a p-type doping effect.^14,15 Recent experiments have shown that boron can be embedded in graphene at higher concentrations (up to 20%).^{10,14,16–18} It was also found that the hexagonal structure exists in boron carbides (B_xC_1−x), thin films with B contents less than 0.5.¹⁹ For B_xC_1−x, thin films with B content at x = 0.25, an ordered layer structure was proposed.^20,21 These make boron-doped graphene (B-graphene) a unique and very attractive material from both fundamental and practical viewpoints. However, the local structural form of B-graphene with relatively high B concentration is not yet well understood. The spatial distribution of doped B in graphene is strongly affected by the synthesis methods and conditions. For the epitaxial B-graphene grown by chemical vapor deposition, the B dopants in graphene can be a preferential substitution of carbon in only one of the graphene sublattices^10,18 and be completely random,²² depending on the metal substrate used. To understand the energetic stability and to find out possible ordered structures of B-graphene, we carry out a global search for the local structure configurations of the B dopants in graphene with a series of concentrations by means of the supercell model and the atomistic simulations.

In the case of substitutional doping of B for C in defect-free graphene, the number of possible structure configurations can be written as $Cn+mm$ if the symmetry is broken, where n and m are the numbers of C and B atoms, respectively. It is evident that the exhaustive search of all possible structure configurations for multiple B dopants in graphene would be prohibitive from the efficiency. To search the stable atomic structures of B-graphene more efficiently, in this work we propose to use a machine learning approach that utilizes Monte Carlo tree search,²³ which was recently utilized for materials and chemical design^24,25 with a rollout schema depending on Bayesian optimization (BO).²⁶ Based on the structure search of B-graphene, we discuss the performance of this proposed scheme of global optimization.

The remainder of this paper is organized as follows. In Sec. II, we briefly introduce the machine learning methods used for the search of atomic structures of B-graphene. The details of computation setup are given in Sec. III. The results for the optimization performance of the machine learning methods, the stabilities of B-graphene, and the electronic structures of B-graphene are presented in Sec. IV. The conclusions are given in Sec. V.

Determination of optimal material structure with certain quality metrics traditionally depended on the experience of domain knowledge experts and trial-and-error experiments. Several machine learning-based methods have been proposed to accelerate this process with as few experiments as possible.^24,27–31 In such methods, the structure prediction is often formulated as a selection of optimal solution from a candidates space that maximizes or minimizes a black-box function (usually the target property).^32,33 Experiments are commonly replaced by simulators such as first-principles calculations.

Currently available efficient methods such as Bayesian optimization (BO)^26,34 are not scalable enough. Due to its exceptional performance in computer Go game,^23,35 Monte Carlo tree search (MCTS)²³ has recently gained attention in material design²⁴ offering a superior scalability with efficiency trade-off.

MCTS employs a shallow search tree where a node is a possible assignment of an atom into a position in the structure. Iteratively, the tree expands towards the optimal solution in 4 steps: selection, expansion, simulation, and backpropagation. A full solution is obtained from this shallow tree using the random rollout (completion) technique in the simulation step. To increase the MCTS efficiency, we propose to use Bayesian learning to engineer the rollout²⁴ instead of the random selection. We briefly discuss this rollout operation here. See the supplementary material for details on the proposed method. The source code is freely available at https://github.com/tsudalab/MDTS.

Bayesian optimization methods maintain a surrogate model of the black-box function, most commonly Gaussian process (GP). A pool of candidate S′ is generated at a node of level p_ℓ where each data point represents a full structure of N positions, with positions p₁, …, p_ℓ being determined along the path from the root to the selected node and positions p_ℓ+1, …, p_N being randomly generated. Data points in S′ are then vectors of binary values (0, 1). If no data point from S′ is previously observed, GP starts with an initial set of randomly selected data points from S′. GP is updated as more data points are observed. An acquisition function is then used to determine the next optimal solution. Within computational budget, the selected solutions are returned (Fig. 1).

To study B-graphene, a supercell constructed by the 4 × 4 extension of the hexagonal unit cell of graphene was employed and the substitution of carbon by boron was considered. To avoid the spurious interaction between graphene layers, a vacuum thickness in the supercell was set to 20.0 Å. The in-plane lattice constant of graphene supercell was fixed at 4a₀, where a₀ = 2.464 Å for the calculated lattice constant of graphene. Figure 2 shows the atomic structure 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 from 3.125% to 31.25%. The number of possible structure configurations for different numbers of doped B atoms in a 4 × 4 graphene supercell is summarized in Table I. Because of the large number of structure configurations for multiple B dopants (n_B $⩾$ 4) in graphene, we employed a combinatorial structure-generation approach for the local-level modeling of atomic substitutions and partial occupancies in crystals to reduce the number of potential configurations.³⁶ Our proposed algorithm is then applied to search the stable structure configurations of B-doped graphene. During the structure search, the atomic positions in each structure configuration were optimized by the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method.³⁷ The energy and atomic force were calculated by the classical force field method,³⁸ as implemented in the ATK-classical simulator of Atomistix ToolKit (ATK).^39,40 The interactions between the atoms in the system under study were described by the Tersoff bond order potentials developed by Matsunaga et al.⁴¹ for C–B and C–C interactions, which have been widely used to study the B-doped carbon nanostructures. For the top ten structures with lower energies found by the global search methods for each case of multiple B dopants, we carried out density functional theory (DFT) calculations to optimize the atomic position and calculate the energy.

Density functional theory calculations have been performed with the Vienna Ab initio Simulation Package (VASP).^42,43 Perdew-Burke-Ernzerh of (PBE) exchange-correlation functional⁴⁴ within the generalized gradient approximations (GGA) was used in conjunction with the projector-augmented wave (PAW) method.^45,46 A plane wave basis set with an energy cutoff of 500 eV was used and a k-point mesh was sampled with a 11 × 11 × 1 Monkhorst-Pack scheme.⁴⁷ The criterions for total energy and force convergence were set to 10⁻⁵ eV and 10⁻² eV/Å, respectively. Because of the fairly delocalized nature for the defect states induced by doped boron atoms,^48,49 the spin polarization of B-graphene is energetically unfavorable, which has been confirmed in the test calculations of spin-polarized DFT. Therefore, we report here the results from the non-spin-polarized calculations. Since the well-known underestimation of band gap in the GGA-PBE calculations due to self-interaction errors, we have also performed the HSE06 hybrid functional calculations⁵⁰ to study the electronic structures of B-graphene in some cases.

As listed in Table I, the number of symmetry non-equivalent atomic structures for multiple doped B atoms considered in a 4 × 4 graphene supercell can reach several hundred thousands. For such large space of candidate structures of B-graphene, the atomic structure search was carried out by the MCTS with Bayesian rollout. We present in Fig. 3 the evolution of lattice energy during the structure search. In the cases considered here, we can see that the change of lattice energies for found optimal structures is less than 0.02 eV per supercell 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 structures, the finding of optimal structure takes 40% of the search space. But for the case of search space with candidate structures more than one hundred thousand, the optimal structure of B-graphene could be found within 2% of the search space. Therefore, our proposed optimization scheme might be very suitable for the larger search space that is hard to be treated, for example, the use of larger supercell size with a kept number of dopant atoms.

To evaluate the stabilities of doped B atoms in graphene, we calculate their formation energies according to the following definition:

where E_t(BG) and E_t(PG) are the total energies of graphene supercell with and without B dopants, m_B is the number of doped B atoms, and μ_B and μ_C are the chemical potentials of boron and carbon, which are taken as the total energies per atom of α-boron crystal and pristine monolayer graphene, respectively. The formation of single B substitution in the 4 × 4 graphene supercell is about 1.186 eV, which is in good agreement with the result (1.12 eV) reported in the previous study.⁵¹ 

It is well known that the honeycomb lattice of graphene consists of two sublattices, namely, A and B sublattices. To understand the occupancy preference of multiple-doped B atoms in graphene, it is necessary to first examine the interaction of two B substitutions in graphene. For this purpose, we calculate the interaction energy of two B substitutions as defined in the following equation:

where E_t(PG), E_t(1B), and E_t(2B) are the total energies of graphene supercells with zero, one, and two doped B atoms, respectively. The positive (negative) sign of the interaction energy defined in Eq. (2) indicates two substitutional B dopants in graphene that repel (attract) each other. It also implies the increase (decreasing) in the formation energy of two B substitutions with respect to that of single B substitution. In the 4 × 4 graphene supercell considered here, the most stable configuration 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. This indicates that the interaction of two substitutional B atoms in graphene is repulsive. For two B substitutions at the first nearest-neighboring (NN) lattice sites (called ortho configuration) in graphene, they show the strongest repulsive interaction and the corresponding interaction energy is 1.262 eV. For two B substitutions at the second and third NN lattice sites (called meta and para configurations, respectively) in graphene, their interaction energies are 0.507 eV and 0.075 eV, respectively. As the para configuration of the B-B pair is more stable than the meta configuration, it has been found in the chemically synthesized B-doped nanographene by using the boron-containing polycyclic aromatic hydrocarbon (PAH).^52,53 The B-graphene with such less stable meta configuration of B-B pair was synthesized by chemical vapor deposition (CVD) on a nickel or cobalt substrate that has strong interaction with graphene.^10,18

The atomic structures for the most stable configurations of B-graphene at different B concentrations are depicted in Fig. 4. They have been confirmed in the calculations that take into account the relaxation of in-plane lattice constants. We have verified their stabilities in the calculations of local-density approximations (LDAs) and meta-GGA. The detailed comparison of the results is given in the supplementary material. The corresponding formation energies of B substitutions are presented in Fig. 5. In these structure configurations, the para configuration of the B-B pair is dominant, indicating the substitutional B atoms prefer to occupy different sublattices in the free standing pristine graphene. The unbalanced sublattice doping observed in experiments may be attributed to the strong interaction of metal substrate and graphene.^10,18 In particular, it is noticed that the formation energy of four B substitutions is lower than those of three and five ones. In the case of eight B substitutions, which corresponds to BC₃, its formation energy is lower than those of seven and nine ones. For the most stable structure configuration of eight B substitutions, as shown in Fig. 4(g), all of the benzene-like C6 rings are separated by B-B pairs and thus its aromaticity is very high. This structure configuration has been proposed in the literature^20,21 for the synthesized BC₃ compound.

To explore the electronic structure of graphene modified by B doping, we present in Fig. 6 the total density of states (DOS) of B-graphenes at a series of B concentrations. It is well known that the Fermi level (E_F) of perfect graphene coincides with the Dirac point. Because boron has one less electron than carbon, the B substitution leads to the shifting of E_F down into the valence bands. We also notice that the B doping opens up a band gap at the Dirac point of graphene, which is supported by the electrical conductivity measurement of B-graphene.¹⁴ In the case of eight B atoms doped in the 4 × 4 graphene supercell, a band gap is open around the E_F and thus this structure of B-graphene is a semiconductor, which is consistent with the previous results.^54,55 The band gaps of BC₃ predicted by the GGA-PBE and HSE06 hybrid functional calculations in the present study are about 0.17 eV and 1.38 eV, respectively. To quantitatively characterize the doping level of carriers, we present the work function of B-graphene in Fig. 7. The work function of perfect monolayer graphene predicted by the GGA-PBE calculations in the present study is about 4.25 eV, which is in good agreement with the experiment results (4.272 eV for CVD-grown graphene measured by UV photoelectron spectroscopy⁵⁶ and 4.57 ± 0.05 eV for mechanically exfoliated graphene measured by a scanning Kelvin probe microscope,⁵⁷ respectively.). It also agrees well with the previous result (4.20 eV) of GGA-PBE calculations in the study of Lazar et al.⁵⁸ We can see that the boron doping increases the work function of graphene significantly. The similar trend was also found in the study of Lazar et al.⁵⁸ Considering the underestimation of band gap in the GGA-PBE calculations, we have also performed the HSE06 hybrid functional calculations to check the work functions of perfect monolayer graphene and B-graphene, which includes the cases of singe B substitution in 4 × 4 supercell and BC₃. The obtained work functions for these three cases are 4.23, 5.23, and 5.80 eV, respectively. Depending on the doped B concentration, the GGA-PBE calculations show that the work function of graphene can be tuned by 0.7 eV–1.0 eV.

In summary, we have employed the Monte Carlo tree search (MCTS) with Bayesian rollout to search the stable structures of B-graphene for the boron concentration up to 31.25%. Compared with the sole use of Bayesian optimization, this integrated optimization method shows a superiority of better scalability. Our results show that in the free-standing graphene, the doped boron atoms energetically prefer to substitute for the carbon atoms at different sublattices. For B-graphene with high boron concentration, the para B-B pairs are dominant. Because of the repulsive interaction between boron substitutions in graphene, the doped boron would tend to be segregated. The doped boron can open a band gap at the Dirac point of graphene. Particularly in the concentration of doped B at 25%, namely, BC₃, the found stable structure exhibits semiconducting behavior. We also find that boron doping can lead to an increase of the work function of graphene in the range between 0.7 eV and 1.0 eV. Our results would be very helpful to further explore the application of B-graphene. In addition, the proposed method can be applied to multiple atom types assignment problems such as boron and nitrogen codoped graphene by allowing more than two branches in the search tree and using a descriptor based on one-hot encoding.

See supplementary material for the details on the proposed machine learning method and the detailed comparison of the results obtained by the different exchange-correlation functionals in DFT calculations.

This work was partially supported by “Materials research by Information Integration” Initiative (MI²I) project of the Support Program for Starting Up Innovation Hub from Japan Science and Technology Agency (JST). The computation in this study was performed on Numerical Materials Simulator at NIMS.

The authors would like to thank Kei Terayama and Shenghong Ju for their fruitful discussion. Also, we are thankful for anonymous referees for their comments and suggestions to improve the manuscript.