Boron adatom adsorption on graphene: A case study in computational chemistry methods for surface interactions

Though weak surface interactions and adsorption can play an important role in plasma processing and materials science, they are not necessarily simple to model. A boron adatom adsorbed on a graphene sheet serves as a case study for how carefully one must select the correct technique from a toolbox of computational chemistry methods. Using a variety of molecular dynamics potentials and density functional theory functionals, we evaluate the adsorption energy, investigate barriers to adsorption and migration, calculate corresponding reaction rates, and show that a surprisingly high level of theory may be necessary to verify that the system is described correctly.

1 Introduction

Examining the adsorption of atomic boron on a graphene sheet offers insight into the difficulties of modelling surface interactions, which can play an important role in materials processing applications. Our interest in using classical molecular dynamics to model the process by which boron is deposited on graphite plasma facing components in fusion devices led us to wonder how boron adheres to the surface of a polycrystalline piece of graphite, though admittedly these adsorption processes occur at far lower energies than those relevant for a fusion device. A wide variety of computational approaches are available to investigate adsorption processes, including classical molecular dynamics (MD), density functional theory (DFT) methods, post-Hartree-Fock methods, and multi-configurational methods. It can be difficult to know ahead of time the right balance between computational cost and the level of accuracy required by the problem under study. Thus, understanding the level of theory necessary to correctly identify binding sites and surface chemistry is important, especially for weak interactions–such as the adsorption of atomic boron on a graphene sheet.

It is worth noting that in addition to general interest as an avenue for investigating methods of modelling surface adsorption, the structure of atomic boron adsorbed on graphene has a variety of suggested applications. It has been proposed as an intermediate structure in a mechanism for the barrier-free substitutional doping of graphene with boron atoms, a material which can act as a sensor for carbon monoxide gas, and a structure with a localized magnetic moment [1–3]. This last application was investigated by Li et al in a study in which the structure was both modeled using DFT calculations and observed experimentally and probed with x-ray photoelectron spectroscopy [3]. Other studies of atomic boron adsorbed onto graphene have used only DFT methods to identify the adsorption site and energy [1, 2, 4–7].

To illustrate the difficulty of achieving some reasonable approximation of chemical accuracy when investigating weak surface interactions such as the adsorption of boron on graphene surfaces, we performed calculations using molecular dynamics (MD) potentials, DFT methods with a variety functionals, and second order Moller-Plesset (MP2) methods.

2 Computational methods

There are three possible binding sites for adatoms on the surface of a graphene sheet: the top site (T-site) directly above a carbon atom, the bridge site (B-site) directly over the bond between two carbon atoms, and the hollow site (H-site) centered over a hexagon. (See Figure 1). Previous DFT studies have found that for a boron adatom, the B-site is the most favorable [1–7]. Calculated adsorption energies and distances between boron and the nearest carbon atom are shown in Table 1.

FIGURE 1

FIGURE 1. Diagram of the three possible adatom binding sites on a graphene sheet.

TABLE 1

TABLE 1. Adsorption energies (E_ads) and nearest boron-carbon distances (d_B-C) for past DFT studies of atomic boron adsorbed on graphene and graphene-like surfaces. Distances marked with an * were inferred assuming a graphene C-C bond distance of 1.42 Å.

Though all previous studies agree that the B-site is the most energetically favorable, the discrepancy in the calculated binding energy is quite striking. The values given range from 0.24 to 1.8 eV, though the distance between the boron atom and nearest carbon atoms of the stable structure, when given, was calculated to be ∼1.8 Å in each of the previous studies [1–7]. Notably, many of these calculations were implemented using the PBE functional of Perdew, Burke and Ernzerhof, with a handful of exceptions [8]. To our knowledge, the only calculation of the energy barrier for the migration of boron adatoms between adjacent binding sites was done by Nakada and Ishii using an LDA functional. They estimated this migration barrier to be 0.12 eV [5].

2.1 Classical molecular dynamics methods

While typically less accurate than quantum chemistry methods, classical molecular dynamics (CMD) allows for the simulation of much larger systems over much larger time scales due to its lower computational cost. Interatomic potentials are used to calculate the forces on particles, and the equations of motion are numerically solved to track the trajectory of the particles over time. For simulations that involve accurately modelling bond breaking and formation, bond order potentials (BOP) enable the description of different bonding states between atoms. These interatomic potentials are produced by parameterizing empirical formulae based on data often produced via quantum chemistry calculations, such as DFT descriptions of bond lengths, bond angles, and the potential energy surface near pertinent molecular geometries. Naively, one might assume that a MD potential optimized on molecules and structures sufficiently similar to the relevant surface interaction would yield a reasonable estimate for the adsorption energy. MD potentials are regularly used to simulate reactions that are not identical to the reactions on which they were fitted. Though this may work well in some situations, for applications involving weak interactions such as the adsorption of boron on graphene, the fidelity of these potentials requires validation.

Our CMD simulations were run using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code [9]. In order to model the adsorption of boron adatoms on graphene, we investigated three different BOPs. These include a Tersoff-type potential and three ReaxFF potentials. [10–13] This particular Tersoff-type potential was developed to investigate thermal conductivity in hexagonal boron nitride and graphene hybrid nanostructures and was parameterized to reproduce DFT energetics obtained using the PBE functional–specifically, the energy as a function of interface distance between hexagonal BN and graphene sheets [10]. The ReaxFF potentials belong to a family of bond order potentials first developed in 2001 and modified to its current form in 2008 [14, 15]. The first ReaxFF potential (ReaxFF-1) was developed to investigate deuterium uptake in amorphous carbon surfaces containing varying amounts of lithium, boron, and carbon, and it is a modified form of a potential initially developed to model chemistry at the anode-electrolyte interface of lithium-sulfur batteries [11, 16]. The data used for parameterization includes DFT energetics of lithium interacting with a boron substituted cyclohexane ring [11]. The second ReaxFF potential (ReaxFF-2) was developed to model shear deformations of boron carbide, and like the Tersoff-type potential, it was also parameterized to match DFT calculations obtained using the PBE functional. The parameters were fitted so that the potential matched data regarding the interaction of two B₁₀C₂H₁₂ icosahedra, the equations of state and heats of formation for various phases of boron and boron carbide, and the shear deformation of boron carbide [12]. The third and final ReaxFF potential (ReaxFF-3) was developed to model liquid CBN (carbon-boron-nitride) hydrogen-storage materials, parameterized by fitting DFT calculations of potential energy surfaces related to varying bond lengths/angles for a variety of small molecules containing carbon, hydrogen, and/or boron [13]. This last ReaxFF potential is in fact part of the “aqueous branch” of the ReaxFF development tree, whereas the first two are from the “combustion branch.” [17] In conjunction with the ReaxFF implementation, the charge equilibration (QEq) method was used to compute the charges on each atom at every timestep [18].

We prepared the initial geometry of the graphene sheet by constructing a bilayer graphite crystal with planar dimensions of approximately 40 Å by 35 Å and periodic boundary conditions. After allowing the system to relax to a minimized energy (including the dimensions of the simulation cell) the top layer of graphene was deleted, leaving a single layer of graphene consisting of 512 atoms. This is the equivalent area of 16 × 16 unit cells of graphene. We used the resulting geometry to investigate graphene-boron surface interactions as the carbon atoms of the graphene sheet were held stationary. For each interatomic potential and each of the three binding sites specified in Figure 1, a single boron atom was moved from 4.5 Å above the binding site towards the graphene sheet and the potential energy of the system was recorded along this trajectory. Minima along each one-dimensional potential energy scan were compared to identify the most stable binding site for each MD potential.

2.2 Quantum chemistry methods

All quantum chemistry calculations were carried out using cluster models in Gaussian 16 software [19]. Large planar aromatic hydrocarbons–bisanthene (C₂₈H₁₄) and coronene (C₂₄H₁₂)—were used as a model for graphene. Such polyaromatic hydrocarbons have previously been used as a model for graphene in computational investigations of adsorption on graphene surfaces [20, 21]. Though adsorption energies calculated using finite models may differ from adsorption energies calculated for an infinite graphite sheet represented by periodic boundary conditions, we expect qualitative agreement regarding the predicted stable geometries. Binding energies calculated for methane interacting with coronene have previously been shown to be within a few tenths of an eV to the binding energies predicted for methane interacting with much larger polyaromatic hydrocarbons (consisting of hundreds of atoms), indicating that coronene is sufficiently large to qualitatively model adsorption on infinite graphene sheets [22].

In addition to adsorption on coronene and bisanthene, we explored the interactions of atomic boron with benzene (C₆H₆). These calculations were performed mainly for comparison with the calculations previously done by Bettinger and Kaiser using multiconfigurational methods to study the formation of benzoborirene [23]. Though benzene is too small to be a good model for an infinite graphene sheet, its interactions with boron serve as a way to compare the performance of various DFT functionals against more accurate calculations performed at a higher level of theory. We used three different DFT functionals of the generalized gradient approximation (GGA) form, including the PBE functional described earlier in this work, Becke’s three-parameter hybrid exchange functional (B3LYP), and a hybrid exchange functional with empirical long-range dispersion (wB97XD) [8, 24, 25]. We employed these functionals in conjunction with the double zeta correlation-consistent Dunning basis set, cc-pVDZ [26]. The default convergence criteria of the Gaussian 16 code were implemented–namely, energy convergence is reached if the energy difference between two self consistent field (SCF) cycles is less than 10⁻⁶ Hartree. Adsorption energies were calculated as below:

${\mathit{E}}_{\mathit{a}\mathit{d}\mathit{s}}={\mathit{E}}_{{\mathit{C}}_{\mathit{x}}{\mathit{H}}_{\mathit{y}}\mathit{B}}-{\mathit{E}}_{{\mathit{C}}_{\mathit{x}}{\mathit{H}}_{\mathit{y}}}-{\mathit{E}}_{\mathit{B}}(1)$

Where $E_{C_x{H}_{y}B}$ represents the electronic energy of the bound structure, $E_{C_x{H}_y}$ represents the electronic energy of the hydrocarbon alone, and $E_B$ the electronic energy of a lone boron atom. Energy barriers are calculated similarly, as the difference between the transition state energy and the energy of the initial configuration, whether that be separate products or the adsorbed state. All energies have been corrected by the zero-point vibrational energy.

Transition states corresponding to barriers to adsorption and barriers to migration between adjacent binding sites were investigated as well. These transition states correspond to saddle points in the potential energy landscape, and are identified by the presence of one imaginary frequency in the vibrational spectrum. Reaction rates were calculated from the associated energy barriers as a function of temperature, T:

${\mathit{k}}_{\mathit{a}\mathit{d}\mathit{s}}=\frac{{\mathit{k}}_{\mathit{B}}\mathit{T}}{\mathit{h}}\frac{{\mathit{Z}}_{\mathit{v}\mathit{i}\mathit{b}}^{\mathit{T}\mathit{S}}}{{\mathit{Z}}_{\mathit{v}\mathit{i}\mathit{b}}^{{\mathit{C}}_{\mathit{x}}{\mathit{H}}_{\mathit{y}}}{\mathit{Z}}_{\mathit{t}\mathit{o}\mathit{t}}^{\mathit{B}}}{\mathit{e}}^{-{\mathit{E}}_{\mathit{b}\mathit{a}\mathit{r}}/{\mathit{k}}_{\mathit{B}}\mathit{T}}(2)$

${\mathit{k}}_{\mathit{m}\mathit{i}\mathit{g}}=\frac{{\mathit{k}}_{\mathit{B}}\mathit{T}}{\mathit{h}}\frac{{\mathit{Z}}_{\mathit{v}\mathit{i}\mathit{b}}^{\mathit{T}\mathit{S}}}{{\mathit{Z}}_{\mathit{v}\mathit{i}\mathit{b}}^{{\mathit{C}}_{\mathit{x}}{\mathit{H}}_{\mathit{y}}\mathit{B}}}{\mathit{e}}^{-{\mathit{E}}_{\mathit{b}\mathit{a}\mathit{r}}/{\mathit{k}}_{\mathit{B}}\mathit{T}}(3)$

Here, ${\mathit{k}}_{\mathit{B}}$ is the Boltzmann constant and $\mathit{h}$ is Planck’s constant. ${\mathit{Z}}_{\mathit{v}\mathit{i}\mathit{b}}^{\mathit{T}\mathit{S}}$, ${\mathit{Z}}_{\mathit{v}\mathit{i}\mathit{b}}^{{\mathit{C}}_{\mathit{x}}{\mathit{H}}_{\mathit{y}}}$, ${\mathit{Z}}_{\mathit{t}\mathit{o}\mathit{t}}^{\mathit{B}}$, and ${\mathit{Z}}_{\mathit{v}\mathit{i}\mathit{b}}^{{\mathit{C}}_{\mathit{x}}{\mathit{H}}_{\mathit{y}}\mathit{B}}$ are the partition functions associated with the vibrational energy of the relevant transition state, the vibrational energy of the hydrocarbon alone, the total energy of the boron atom, and the vibrational energy of the bound state. E_bar is the height of the energy barrier for the relevant process.

3 Results and discussion

3.1 Classical molecular dynamics results and discussion

The potential energy curves are presented in Figure 2, scanning over the distance between the boron atom and the graphene sheet. For each MD potential, the preferred binding site and corresponding binding energy minimum are presented in Table 2. The adsorption energies predicted by the ReaxFF potentials are larger than any of those predicted by DFT methods in Table 1, and the adsorption energy of the Tersoff-type potential is lower. Examining the potential energy curves for the ReaxFF potentials we find large and unphysical fluctuations in the potential energy as the boron atom moves closer to the graphene sheet. The Tersoff-type potential and two of the ReaxFF potentials predict that the hollow site is the most stable, whereas the other ReaxFF potential predicts that the bridge site is the most energetically favorable position for the boron adatom.

FIGURE 2

FIGURE 2. Potential energy as a function of boron graphene distance (z_B/G) for the (A) hollow (B) top, and (C) bridge sites. The energy of the Tersoff-type potential [10] is scaled by a factor of 10 to improve visibility of the minimum in comparison to those of the ReaxFF potentials, ReaxFF-1 [11, 16], ReaxFF-2 [12], and ReaxFF-3 [13].

TABLE 2

TABLE 2. Adsorption energy and preferred binding site for investigated MD potentials.

None of these MD potentials were parameterized to reproduce the energetics of a boron adatom adsorbing on a graphene sheet, and their unsuitability for this task is not a mark against them. The systems used to parameterize the potentials often did not remotely resemble graphene, with the exception of the Tersoff-type potential. However, we can see that even the Tersoff-type potential predicted a binding energy that fell outside the broad range of binding energies predicted by previous quantum chemistry calculations. The two potentials that were optimized based on DFT calculations for smaller molecules (ReaxFF-1 and ReaxFF-3) predicted unreasonably large adsorption energies, whereas the potentials that were optimized based on DFT calculations for systems containing slabs of graphene and boron carbide (Tersoff-type and ReaxFF-2, respectively) were smaller and more reasonable. Still, the potential energy scans sufficiently demonstrate that the investigated carbon-boron bond order potentials are unsuitable for reproducing the surface interactions of boron adatoms adsorbing on graphene, and one should not assume that molecular dynamics potentials will be capable of accurately modelling surface interactions if they have not been optimized using DFT calculations specific to those reactions.

3.2 Quantum chemistry results and discussion

3.2.1 Large planar hydrocarbons interacting with boron

It is immediately apparent from the optimized geometries shown in Figure 3 that the PBE functional predicts a different binding site than the other two functionals. While B3LYP and wB97XD predict that the boron atom will bind near a top site, canted slightly towards a carbon-carbon bond, the PBE functional predicts that the boron atom will bind at the bridge position, as it has in previous studies. The predicted adsorption energies for boron on bisanthene using the B3LYP, wB97XD, and PBE functionals were −0.72 eV, −1.02 eV, and −1.11 eV, respectively. (See Table 3). These fall within the range of previously predicted adsorption energies of boron on graphene. (See Table 1).

FIGURE 3

FIGURE 3. Binding sites on C₂₈H₁₄ calculated using (A) B3LYP/cc-pVDZ (B) wB97XD/cc-pVDZ, and (C) PBE/cc-pVDZ methods. Carbon = Grey, Hydrogen = Black, and Boron = Blue.

TABLE 3

TABLE 3. Adsorption energies (E_ads), barriers to adsorption (E_bar), and barriers to migration (E_mig) calculated with and without zero-point energy corrections for each of the DFT functionals investigated.

Transition states representing barriers to adsorption were identified in for the B3LYP and wB97XD functionals, in addition to a weakly physisorbed state at boron-carbon distances beyond the transition state. Potential energy scans following the intrinsic reaction coordinate path from the weakly physisorbed state, over the transition state, and to the fully adsorbed state are depicted in Figure 4. When zero-point energy (ZPE) corrections are taken into account, the transition state energy predicted with the B3LYP method actually lies below the energy of the separated products. It does not, therefore, pose a barrier to adsorption. In contrast the ZPE-corrected barrier predicted by wB97XD is 0.16 eV. Calculations using the PBE functional located no equivalent transition state, as the potential energy was found to monotonically increase as the boron atom was pulled farther from the surface. Thus, only the wB97XD functional predicts the existence of a barrier to adsorption.

FIGURE 4

FIGURE 4. IRC scans following the reaction from the weakly physisorbed state (left) over the transition state (center) to the fully adsorbed state (right) calculated with the (A) B3LYP/cc-pVDZ method and (B) wB97XD/cc-pVDZ method. The energies depicted are in reference to the energy of the fully separated products, and no zero-point energy corrections have been made. Carbon = Grey, Hydrogen = Black, and Boron = Blue.

Diffusion along the graphene sheet was also considered. In migrating between adjacent binding sites, the boron adatom experiences a barrier to diffusion. This migration was modeled using bisanthene for the B3LYP and wB97XD functionals and coronene for the PBE functional. The predicted barriers to diffusion are 0.20 eV for wB97XD, 0.11 eV for B3LYP, and 0.09 eV for PBE. The potential energy scans along the intrinsic reaction coordinates for this migration process can be seen in Figure 5.

FIGURE 5

FIGURE 5. IRC potential energy scans for the migration of the boron adatom between adjacent binding sites, using the (A) B3LYP functional (B) wB97XD functional, and (C) PBE functional. Scans performed in (A) and (B) used a bisanthene molecule as a model for graphene, the scan performed in (C) used coronene. The energies depicted are in reference to the energy of the fully bound state, and no zero-point energy corrections have been made. Carbon = Grey, Hydrogen = Black, and Boron = Blue.

The reaction rates calculated from the variety of energy barriers described previously are plotted as a function of temperature in Figure 6. A summary of the various calculated energies (adsorption energies, barriers to adsorption, and barriers to migration) with and without zero-point energy corrections can be found in Table 3.

FIGURE 6

FIGURE 6. Reaction rates as a function of temperature for the processes of (A) boron adatom adsorption on graphene and (B) migration of the boron adatom between adjacent binding sites on the graphene surface.

3.2.2 Benzene interacting with boron

Given the discrepancies found above, it seems that we must inevitably conclude that either the PBE functional or both the B3LYP and wB97XD functional incorrectly predict the binding site for boron adsorption on graphene. One would tentatively expect the hybrid functionals to be more accurate than the PBE functional in predicting reaction barriers and binding energies, existing on a “higher rung” of the DFT functional “Jacob’s ladder” to heavenly chemical accuracy, the popular analogy first used by Perdew [27]. However, for further assurance we turn to modeling the interaction of atomic boron with a smaller aromatic planar hydrocarbon–benzene. Though the interaction between benzene and atomic boron is a crude model for the interaction between graphene and atomic boron, the six-membered carbon ring of benzene is nevertheless a sufficient proxy for the six-membered ring of graphene. Both benzene and graphene are aromatic compounds with delocalized electrons shared among pi bonds; the main difference is that the conjugated system of a benzene ring is terminated by hydrogen atoms, whereas in graphene the conjugated system of a six-membered ring is part of the graphene sheet. Therefore, we can use this very small system to compare the results of different DFT functionals to more accurate (and more computationally expensive) calculations done at a higher level of theory.

This comparison will enable us to make reasonable predictions about which DFT functionals are producing accurate results for the interactions of large polyaromatic hydrocarbons with atomic boron.

We once again see that the PBE functional predicts a stable structure with the boron atom positioned directly over the bond between two carbon atoms, whereas the B3LYP and wB97XD functionals both predict a stable geometry with the boron atom nearer to one of the carbon atoms. (See Figure 7). These geometries can be directly compared to structures described by Bettinger and Kaiser, who used B3LYP and a complete active space (CAS) method to investigate the formation of benzoborirene. [20]. Despite using a different basis set, the stable geometry predicted by their B3LYP method agrees with the stable geometry predicted by our B3LYP method. More importantly, the CAS calculations of Bettinger and Kaiser indicate that a stable structure exists when the boron is directly over one of the carbon atoms (analogous to T-site adsorption on graphene), and that the geometry predicted to be stable by the PBE functional (labeled 5 by Bettinger and Kaiser) is actually a transition state.

FIGURE 7

FIGURE 7. Stable geometries of atomic boron interacting with benzene as predicted by (A) B3LYP/cc-pVDZ calculations and (B) PBE/cc-pVDZ calculations. Carbon = Grey, Hydrogen = White, and Boron = Pink.

As the multi-configurational complete active space self-consistent field (CASSCF) method used in this paper is at a much higher level of theory than any of our density functional theory methods, this indicates that the B3LYP and wB97XD hybrid functionals predict the binding geometry with greater accuracy than the PBE functional. Since many studies of adatom adsorption on graphene have made use of the PBE functional, further investigation of these structures using hybrid functionals such as B3LYP and wB97XD may produce better estimates of adsorption energy and possibly predict different binding sites than those predicted in previous literature.

Interestingly, MP2 calculations did not agree with the CASSCF results, as this post-Hartree-Fock method agreed with the PBE functional in its prediction that the structure shown in Figure 7B was stable. Typically, the MP2 method is used to validate DFT results [28]. However, it has been shown elsewhere that the interaction of a beryllium atom with a benzene ring involves multi-reference electronic states which MP2 fails to account for, and the MP2 method overestimates the adsorption energy of this system [29]. We might reasonably anticipate that the MP2 method also overestimates the attractive interactions between a boron atom and a benzene ring, and we would generally expect MP2 to be less accurate than the more sophisticated CAS methods used by Bettinger and Kaiser [23].

4 Conclusion

We have investigated the adsorption of atomic boron on graphene, as modeled by bond order MD potentials, and as modeled by large planar hydrocarbons using a variety of DFT functionals. The MD potentials considered here are not suitable for modeling the adsorption of a boron adatom on a graphene sheet, as one might expect considering that they were not parameterized with this interaction in mind. However, the very unphysical oscillations in some of the ReaxFF potentials and the magnitude of the binding energies predicted should serve as a reminder that MD potentials used to investigate processes where weak surface interactions play a role ought to be optimized to reproduce the energetics of the relevant reactions. To our knowledge, there is no such potential for modeling boron adatom adsorption on graphene. In the future, MD potentials that harness machine learning via the SNAP and DeePMD formalisms may be worth exploring for modelling such surface interactions [30, 31].

Three different DFT functionals were used to estimate binding energies, barriers to adsorption, and barriers to migration between adjacent binding sites for boron adatoms on graphene. For the predicted energy barriers, corresponding reaction rates were calculated. We suggest that the bridge binding site predicted by previous studies may be incorrect, and that the true binding site may be the top site. If this is the case, then the PBE functional may yield the wrong adsorption site, and this is important to note not only for its own sake but for the future parameterization of MD potentials dealing with adsorption on surfaces, which rely on DFT calculations to provide a sufficiently realistic description of relevant potential energy surfaces. Care must be taken to correctly assess the binding energies and barriers related to weak surface interactions, and a surprisingly high level of theory may be needed to verify the correct binding site.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

All authors contributed to the conception of this study. SJ ran the quantum chemistry calculations and wrote the first full draft of the manuscript. AR ran the classical molecular dynamics calculations and wrote the first draft of the corresponding sections of the manuscript. YB provided insight on quantum chemistry calculations. SE helped with installation of software and compiling the codes. IK organized the research effort and managed group interactions. All authors contributed to manuscript revision and approved the submitted version.

Funding

The research described in this paper was conducted under the Laboratory Directed Research and Development (LDRD) Program at Princeton Plasma Physics Laboratory, a national laboratory operated by Princeton University for the U.S. Department of Energy under Prime Contract No. DE-AC02-09CH11466. The United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. In addition, this research used computing resources on the Princeton University Adroit Cluster and PPPL cluster.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Pontes R, Fazzio A, Dalpian G. Barrier-free substitutional doping of graphene sheets with boron atoms: Ab initio calculations. Phys Rev B (2009) 79(3):033412. doi:10.1103/PhysRevB.79.033412

CrossRef Full Text | Google Scholar

2. Fani S, Tavangar Z, Kazempour A. Boron decorated graphene nanosheet as an ultrasensitive sensor: the role of coverage. J Mol Model (2019) 25(6):166. doi:10.1007/s00894-019-4046-z

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Li Q, Lin H, Lv R, Terrones M, Chi L, Hofer W, et al. Locally induced spin states on graphene by chemical attachment of boron atoms. Nano Lett (2018) 18(9):5482–7. doi:10.1021/acs.nanolett.8b01798

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Kotakoski J, Krasheninnikov A, Ma Y, Foster A, Nordlund K, Nieminen R, et al. B and N ion implantation into carbon nanotubes: Insight from atomistic simulations. Phys Rev B (2005) 71(20):205408. doi:10.1103/PhysRevB.71.205408

CrossRef Full Text | Google Scholar

5. Nakada K, Ishii A. Migration of adatom adsorption on graphene using DFT calculation. Solid State Commun (2011) 151(1):13–6. doi:10.1016/j.ssc.2010.10.036

CrossRef Full Text | Google Scholar

6. Pašti I, Jovanović A, Dobrota A, Mentus S, Johansson B, Skorodumova N, et al. Atomic adsorption on pristine graphene along the Periodic Table of Elements – from PBE to non-local functionals. Appl Surf Sci (2018) 436:433–40. doi:10.1016/j.apsusc.2017.12.046

CrossRef Full Text | Google Scholar

7. Yan W, Li X, Zhang X, Cao X, Deng M. Family-dependent magnetism in atomic boron adsorbed armchair graphene nanoribbons. J Mater Chem C Mater (2019) 7(21):6241–5. doi:10.1039/c9tc00186g

CrossRef Full Text | Google Scholar

8. Perdew J, Burke K, Ernzerhof M. Generalized gradient approximation made simple. Phys Rev Lett (1996) 78(18):3865–8. doi:10.1103/PhysRevLett.77.3865

CrossRef Full Text | Google Scholar

9. Plimpton SFast parallel algorithms for short-range molecular dynamics. J Comput Phys (1995) 117(1):1–19. doi:10.1006/jcph.1995.1039

CrossRef Full Text | Google Scholar

10. Kınacı A, Haskins JB, Sevik C, Çağın T. Thermal conductivity of BN-C nanostructures. Phys Rev B (2012) 86(11):115410. doi:10.1103/physrevb.86.115410

CrossRef Full Text | Google Scholar

11. Domínguez-Gutiérrez F, Krstić P, Allain J, Bedoya F, Islam M, Lotfi R, et al. Deuterium uptake and sputtering of simultaneous lithiated, boronized, and oxidized carbon surfaces irradiated by low-energy deuterium. J Appl Phys (2018) 123(19):195901. doi:10.1063/1.5026415

CrossRef Full Text | Google Scholar

12. An Q, Goddard WA III. Atomistic origin of brittle failure of boron carbide from large-scale reactive dynamics simulations: Suggestions toward improved ductility. Phys Rev Lett (2015) 115(10):105501. doi:10.1103/PhysRevLett.115.105501

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Pai SJ, Yeo BC, Han SS. Development of the ReaxFFCBN reactive force field for the improved design of liquid CBN hydrogen storage materials. Phys Chem Chem Phys (2016) 18(3):1818–27. doi:10.1039/C5CP05486A

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Van Duin A, Dasgupta S, Lorant F, Goddard W. ReaxFF:  A reactive force field for hydrocarbons. J Phys Chem A (2001) 105(41):9396–409. doi:10.1021/jp004368u

CrossRef Full Text | Google Scholar

15. Chenoweth K, Van Duin A, Goddard W. ReaxFF reactive force field for molecular dynamics simulations of hydrocarbon oxidation. J Phys Chem A (2008) 112(5):1040–53. doi:10.1021/jp709896w

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Islam M, Bryantsev VS, van Duin ACT. ReaxFF reactive force field simulations on the influence of teflon on electrolyte decomposition during Li/SWCNT anode discharge in lithium-sulfur batteries. J Electrochem Soc (2014) 161(8):E3009–14. doi:10.1149/2.005408jes

CrossRef Full Text | Google Scholar

17. Senftle T, Hong S, Islam M, Kylasa S, Zheng Y, Shin Y, et al. The ReaxFF reactive force-field: Development, applications and future directions. Npj Comput Mater (2016) 2:15011. doi:10.1038/npjcompumats.2015.11

CrossRef Full Text | Google Scholar

18. Aktulga HM, Fogarty JC, Pandit SA, Grama AY. Parallel reactive molecular dynamics: Numerical methods and algorithmic techniques. Parallel Comput (2012) 38(4):245–59. doi:10.1016/j.parco.2011.08.005

CrossRef Full Text | Google Scholar

19. Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, et al. Gaussian 16, Revision A.03. Wallingford CT: Gaussian, Inc. (2016).

Google Scholar

20. Hutama A, Hijikata Y, Irle S. Coupled cluster and density functional studies of atomic fluorine chemisorption on coronene as model systems for graphene fluorination. J Phys Chem C (2017) 121(27):14888–98. doi:10.1021/acs.jpcc.7b03627

CrossRef Full Text | Google Scholar

21. Varenius E. Hydrogen adsorption on graphene and coronene: A van der Waals density functional study. master’s thesis. Göteborg, Sweden: Chalmer’s University of Technology (2011). p. 45.

22. Lazare J, Daggag D, Dinadayalane T. DFT study on binding of single and double methane with aromatic hydrocarbons and graphene: stabilizing CH…HC interactions between two methane molecules. Struct Chem (2021) 32:591–605. doi:10.1007/s11224-020-01657-y

CrossRef Full Text | Google Scholar

23. Bettinger H, Kaiser R. Reaction of benzene and boron atom: Mechanism of formation of benzoborirene and hydrogen atom. J Phys Chem A (2004) 108(21):4576–86. doi:10.1021/jp0375259

CrossRef Full Text | Google Scholar

24. Becke AD. Density-functional thermochemistry. III. The role of exact exchange. J Chem Phys (1993) 98(7):5648–52. doi:10.1063/1.464913

CrossRef Full Text | Google Scholar

25. Chai JD, Head-Gordon M. Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections. Phys Chem Chem Phys (2008) 10(44):6615. doi:10.1039/B810189B

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Dunning TH. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J Chem Phys (1989) 90(2):1007–23. doi:10.1063/1.456153

CrossRef Full Text | Google Scholar

27. Perdew JP, Schmidt K. Jacob’s ladder of density functional approximations for the exchange-correlation energy. AIP Conf Proc (2001) 577:1–20. doi:10.1063/1.1390175

CrossRef Full Text | Google Scholar

28. Porsev V, Barsukov Y, Tulub A. Systematic quantum chemical research of the reaction of magnesium clusters with organic halides. Comput Theor Chem (2012) 995:55–65. doi:10.1016/j.comptc.2012.06.030

CrossRef Full Text | Google Scholar

29. Fernandez N, Ferro Y, Carissan Y, Marchois J, Allouche A. The interaction of beryllium with benzene and graphene: A comparative investigation based on DFT, MP2, CCSD(T), CAS-SCF and CAS-PT2. Phys Chem Chem Phys (2014) 16(5):1957–66. doi:10.1039/c3cp54062f

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Wood MA, Thomspon AP. Extending the accuracy of the SNAP interatomic potential form. J Chem Phys (2018) 148:241721. doi:10.1063/1.5017641

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Lu D, Wang H, Chen M, Lin L, Car R, Weinan E, et al. 86 PFLOPS Deep Potential Molecular Dynamics simulation of 100 million atoms with ab initio accuracy. Comput Phys Commun (2021) 259:107624. doi:10.1016/j.cpc.2020.107624

CrossRef Full Text | Google Scholar