Fe@B6H6 aggregates: from simple building blocks to graphene analogue

We suggest the possibility to build graphene analogue with the planar hexacoordinate wheel-type Fe@B6H6 cluster as the building block through studying theoretically the geometry, stability, and electron structure of its dimer and trimer as well as the dimerization of the two trimers. Employing the dehydrogenation route to polymerization, we can obtain the hexagonal boron sheet that are partly and uniformly filled by Fe atoms in the center of the holes, achieving uniform chemical doping and a very large hexagonal-hole density. Thus, we may offer a novel cluster-assembled material for experimental chemists to construct graphene analogue.


Introduction
Graphene, the first perfect monatomic two-dimensional carbon crystal, was isolated successfully from graphite in 2004 [1]. Its excellent properties, such as extremely high carrier mobility, high thermal conductivity, and high specific surface area, have sparked not only the intensive studies on its synthesis and functionalized applications [2, 3] but also the extensive discoveries towards graphene analogous [4,5], consisting of compositions other than carbon. Boron clusters have been found theoretically [6-9] and experimentally [10-12] to possess a two-dimensional (quasi)planar structure termed boron sheet. However, boron cannot form the stable honeycomb hexagonal-hole framework as graphene because of its electron-deficient character. Instead, part of the hexagonal holes need to be filled by boron atoms adopting buckled form so that the boron sheet can maintain its (quasi)planar structure [11,12]. Therefore, the boron atoms filling the holes serve as an important role in boron sheet and obviously distinguish from the ones forming the hole framework. Under the context, can we choose other atom to partially fill the holes and achieve the better results than that from boron atom? If possible, the more crucial question is how we can fill the holes in boron sheet. The challenging and appealing questions arouse our research interest.
The (quasi)planar hexacoordinate wheel-type clusters, such as CB 6 [2-13] and X@B 6 H 6 (X = V, Cr, and Mn) [14,15], give us some inspiration to answer the questions. The existence of these clusters indicates that it is possible to fill the hole with atoms other than boron. And then, we can use these clusters as building blocks for constructing graphene analogue through polymerization. As we have shown in our recent work [16], the planar wheel-type D 6h Fe@B 6 H 6 with good chemical stability is the global minimum isomer and is therefore more attractive candidate for cluster-assembled materials. Here, we studied theoretically the geometry, stability, and electron structure of the dimer and trimer of Fe@ B 6 H 6 as well as the graphene analogue FeB 6 . Our calculated results indicate that the structures and charges of the dimer and trimer of Fe@B 6 H 6 are similar to those of the monomer Fe@B 6 H 6 . Based on π-electron molecular orbital, Huckel rule, and nucleus-independent chemical shifts, the trimer of Fe@B 6 H 6 can be considered the triphenylene analogue.
Moreover, the phonon dispersion indicates that the graphene analogue FeB 6 has good dynamical stability. These results indicate that the Fe@B 6 H 6 can be used as the building block to build the graphene analogue FeB 6 .

Computational methods
The structures (see Fig. 1) of the monomer, dimer, and trimer of Fe@B 6 H 6 were optimized by employing B3LYP [17,18] method, as implemented in Gaussian 09 [19], with 6-311 + + G(d,p) basis sets for all atoms. To justify the B3LYP method, BP86 [20]/6-311 + + G(d,p) and TPSSh [21]/6-311 + + G(d,p) methods were also used to optimize these structures since the two functionals are regarded as accurate methods for computing structures of transitionmetal complexes [22][23][24]. No imaginary frequency was found in all the calculations. The optimized geometric parameters obtained from the three methods are listed in Table S1. It can be seen that the bonds calculated by the three methods are almost the same. Moreover, most bonds obtained from B3LYP are between the bonds obtained from BP86 and bonds obtained from TPSSh, which indicates that the B3LYP method is reliable for our systems. Therefore, unless otherwise specified, the following studies are carried out by B3LYP method. Using the optimized geometries, the HOMO-LUMO energy gap (∆ H-L ), vertical ionization potential (VIP) [16], vertical electron affinity (VEA) [16], natural population analysis (NPA) [25], and nucleus-independent chemical shifts (NICS) [26] were carried out at B3LYP/6-311 + G(d,p) level. The VIP of monomer, dimer, and trimer is defined by the energy difference between the cationic E(monomer/dimer/trimer) + and neutral E(monomer/ dimer/trimer) calculated at the equilibrium neutral geometry, and the VEA of monomer, dimer, and trimer is the energy difference between the neutral E(monomer/dimer/trimer) and the anionic E(monomer/dimer/trimer) − calculated at the equilibrium neutral geometry, as given in the following equations: For graphene analogue FeB 6 , its optimized structure is carried out by the DMol 3 module implemented in Material Studio 2018 using Perdew-Burke-Ernzerhof (PBE) [27] generalized gradient approximation (GGA) and double-numerical properties plus polarization (DNP). In the convergence Fig. 1 The optimized geometries of D 6h monomer Fe@ B 6 H 6 , D 2d dimer (Fe@B 6 ) 2 H 10 , and C 2v trimer (Fe@B 6 ) 3 H 12 tolerance, the energy, force, and displacement were set as 10 −5 Ha, 0.002Ha/Å, and 0.005 Å, respectively [28]. A vacuum layer of 20 Å is added to avoid the influence of periodic adjacent layers. And Monkhorst-Pack k-mesh of 6 × 6 × 1 is adopted in Brillouin zone.

Results and discussion
Optimized structures and some parameters of the monomer, dimer, and trimer of Fe@B 6 H 6 Figure 1 shows the optimized geometries of Fe@B 6 H 6 and the dimer (Fe@B 6 ) 2 H 10 and the trimer (Fe@B 6 ) 3 H 12 . The bond lengths and the charges of monomer, dimer, and trimer are listed in Tables 1 and 2, respectively. The monomer Fe@ B 6 H 6 forms the perfect regular hexagon with the planar hexacoordinate Fe atom, possessing the highest D 6h symmetry. And the bond lengths and charges agree well with those obtained at BP86/6-311 + G(3df,3pd) level [16]. For the dimer (Fe@B 6 ) 2 H 10 , it has the highest D 2d symmetry and the dihedral angle between the two monomers is 90° and the distance between them is 1.654 Å which is in good agreement with the experimental value of B-B bond length (1.691 Å) [29], suggesting the interaction between the monomers is very strong. While comparing with the monomer Fe@B 6 H 6 , the bond lengths and charges for the dimer (Fe@B 6 ) 2 H 10 do not change significantly in each monomer except the B atom that links the other monomer. For the dimer (Fe@B 6 ) 2 H 10 , all atoms in each monomer are coplanar, which indicates that the character of monomer Fe@B 6 H 6 is well maintained during the dimerization. The dimer (Fe@B 6 ) 2 H 10 with D 2h symmetry is a transition state for the conversion to the D 2d conformation. The energy of D 2d conformation is 3.08 kcal/ mol lower than that of D 2h conformation. To find the barrier separating equivalent D 2d conformers, potential energy curves (see Figure S1) for internal rotation of the dimer were obtained as a function of dihedral angle from − 90 to 90°, in steps of 10°. It can be seen from Figure S1 that the difference energy between the D 2d conformation and the D 2h conformation is the rotation barrier (3.08 kcal/mol). For the trimer (Fe@B 6 ) 3 H 12 , it has the highest C 2v symmetry and the three monomers are perfectly coplanar. The hole in the trimer is not regular hexagon and it is composed of two different types of B-B distances, which are R B1-B5 = 1.665 Å and R B5-B6 = 1.854 Å, respectively. Besides, we also consider the case of non-planar trimers ((Fe@B 6 ) 3 H 14 ). The optimized structures and the relative energies of (Fe@B 6 ) 3 H 14 are shown in Figure S2. To estimate the stability of (Fe@ B 6 ) 3 H 12 , the energy change of the following process (Fe@ B 6 ) 3 H 12 + H 2 → (Fe@B 6 ) 3 H 14 (corresponding to 3 in Figure S2) is calculated. We find that the reaction is endothermic (∆E = 7.81 kcal/mol), indicating that the (Fe@B 6 ) 3 H 12 is thermodynamically stable.
The NPA atomic charges of the monomer, dimer, and trimer of Fe@B 6 H 6 are listed in Table 2. For Fe@B 6 H 6 , the charges of Fe and B are similar to Ref. 16. For (Fe@B 6 ) 2 H 10 and (Fe@B 6 ) 3 H 12 , the charges of Fe, B2, B3, and B4 are almost same as those of Fe@B 6 H 6 . While the charges of B1, B5, and B6 are slightly less than those of Fe@B 6 H 6 , which is attributed to the B-B bond formation. The ∆ H-L , VIP, and VEA of the monomer, dimer, and trimer of Fe@B 6 H 6 are also listed in Table 2. It can be seen that the ∆ H-L values of monomer, dimer, and trimer are 3.58, 3.22, and 2.80 eV, respectively. Although the ∆ H-L value of trimer is the smallest among them, it is larger than the ∆ H-L value (2.63 eV) of triphenylene [30], indicating the monomer, dimer, and trimer of Fe@B 6 H 6 are chemically stable. The VIP values of monomer, dimer, and trimer are 8.06 eV, 7.68 eV, and 7.57 eV, respectively while the VEA values of monomer, dimer, and trimer are 1.69 eV, 2.15 eV, and 2.53 eV, respectively.
For (Fe@B 6 ) 2 H 10 and (Fe@B 6 ) 3 H 12 , the interaction energies (∆E) of monomers are important for their stability. Considering basis set superposition error (BSSE), the ∆E Table 1 The bond distances (Å) in the monomer, dimer, and trimer of Fe@B 6 H 6 at B3LYP/6-311 + + G(d,P)  Table 2 The NPA atomic charges (q/|e|), HOMO-LUMO energy gap (∆ H-L /eV), vertical ionization potential (VIP/eV), and vertical electron affinity (VEA/eV) of the monomer, dimer, and trimer of Fe@B 6 H 6 The molecular orbital and aromaticity of the trimer (Fe@B 6 ) 3 H 12 Since Fe@B 6 H 6 exhibits the similar π molecule orbitals to benzene [16], the trimer (Fe@B 6 ) 3 H 12 may be the triphenylene analogue. In order to confirm our conjecture, the π-electron molecular orbitals (MOs) of (Fe@B 6 ) 3 H 12 and triphenylene are plotted in Fig. 2. It can be seen that the shape of these Mos of (Fe@B 6 ) 3 H 12 and triphenylene is similar. For example, the HOMO-5 of (Fe@B 6 ) 3 H 12 is bond MO which is similar to the bond MO of HOMO-9 in triphenylene. In addition, both (Fe@B 6 ) 3 H 12 and triphenylene have three degenerate MOs (HOMO, HOMO-6, and HOMO-8 in (Fe@B 6 ) 3 H 12 and HOMO, HOMO-2, and HOMO-6 in triphenylene). As a result, their nine π-electron MOs accommodate 18 π electrons that satisfy the (4n + 2) Huckel rule. Thus, the trimer (Fe@B 6 ) 3 H 12 exhibits the aromaticity and can be considered to be the triphenylene analogue. NICS is a simple and efficient criterion to characterize aromatic nature. To better understand the aromaticity, the calculated NICS(d) (d = 0 and 1 for inside and above the hole, respectively.) of (Fe@B 6 ) 3 H 12 and triphenylene are also shown in Fig. 2. The NICS(0) = − 0.53 ppm and NICS(1) = − 0.28 ppm of the hole in the trimer (Fe@ B 6 ) 3 H 12 are less negative than the NICS(0) = − 1.72 ppm and NICS(1) = − 5.09 ppm of the hole in triphenylene, which indicates that the hole of trimer (Fe@B 6 ) 3 H 12 is less aromatic than that of triphenylene. While the monomer in the trimer (Fe@B 6 ) 3 H 12 has very strong aromatic character since its NICS(1) = − 15.2 ppm is more negative than − 9.8 ppm for the monomer in triphenylene, which can compensate the aromaticity of the hole in trimer (Fe@B 6 ) 3 H 12 .

The geometric and electronic structures and stability of graphene analogue FeB 6
Before building the graphene analogue FeB 6 , we examined the bigger stable aggregates. Two kinds of different dimerization of the trimer are shown in Fig. 3. The six monomers reveal perfect coplanarity in each of them, indicating the trimer possesses very good ability of plane expansion. Thus, assembling the stable trimers (Fe@B 6 ) 3 H 12 can provide the possibility of building graphene analogue FeB 6 as the triphenylene in graphene [34,35].
And then, we optimized the graphene analogue FeB 6 , as shown in Fig. 4A. The FeB 6 with D 6h is a completely planar structure which is the third most stable structure of the monolayer FeB 6 [36]. The boron ring with Fe atom in the FeB 6 has the B-B bond length of 1.860 Å and B-Fe bond length of 1.860 Å, which are slightly longer than those (1.827 Å) of Fe@B 6 H 6 monomer. And the bond lengths of two different B-B in the boron ring without Fe atom are 1.860 Å and 1.661 Å, which are similar to those of boron ring in Fe@ B 6 H 6 trimer. The bond lengths in parenthesis are obtained from PW91 [37]/DNP approach. It can be seen that both PBE and PW91 can give the similar results, which indicates that the PBE approach is reliable. Therefore, the graphene analogue FeB 6 preserves the structural features of monomer and trimer of Fe@B 6 H 6 . The electron localization function (ELF) of the FeB 6 is shown in Fig. 4B. It can be seen that the electrons are localized around the B-B bonds; especially, the B-B bonds between Fe@B 6 present highly localized electron distributions, while the center of B-ring and Fe atom exhibits highly delocalized electron distributions.
We also studied the hexagon hole density of the FeB 6 . The hexagon hole density (η) is defined as [30,38]: According to the formula, the triangular boron sheet has η = 0, the hexagonal boron sheet η = 1/3 [38]. For the FeB 6 , it represents a hexagonal-hole density of η = 2/7, which is bigger than those in pure boron α and β [38] and very close to the hexagonal boron sheet η = 1/3.
We also calculated the band structure and the total/ partial density of states (T/PDOS) of FeB 6 monolayer, as shown in Fig. 5. The FeB 6 monolayer can be regarded as a narrow band gap semiconductor because it has a small band gap of 0.4 eV. The analysis of PDOS reveals that the top of valence bands mainly comes from the contributions of Fe-3d orbital and B-2 s orbital and the bottom of No. of hexagon holes No. of atoms in the original triangular sheet conduction bands mainly comes from the contributions of Fe-3d orbital and B-2p orbital. Besides, we also investigated the dynamical stability of the FeB 6 . The phonon dispersion is shown in Fig. 6. The unit cell of FeB 6 monolayer has seven atoms, suggesting that the phonon band structures should have 21 phonon branches. The highest frequency reaches up to 1204 cm −1 , and is higher than the highest frequency of 1036 cm −1 in BSi [39] and 924 cm −1 in Ti 2 B 2 [40], indicative of robust Fe-B and B-B interactions in FeB 6 monolayer. Furthermore, the absence of virtual frequencies at any high-symmetry direction also confirms the dynamic stability of the FeB 6 .
To investigate the thermal stability of FeB 6 monolayer, the NVT molecular dynamics (MD) simulations for 1 × 1 FeB 6 unit cell and 2 × 2 FeB 6 super cell were performed at 5 ps interval at 300 and 600 K with the Nose-Hoover thermostar [41]. The results of MD simulations at 300 K and 600 K are shown in Figure S3. It can be seen that no large fluctuations of the potential energy are found in MD simulations, and the structures of FeB 6 monolayer almost keep the original geometries without any broken bonds. These results confirm that the FeB 6 monolayer is thermally stable. To explore the possibility of the graphene analogue FeB 6 via mechanical stripping, we calculated its exfoliation energy. Here, the n-layer exfoliation energy per unit area is determined by the following formula [42] where E iso (n) and E bulk are the unit cell energy of isolated n-layer slab in vacuum and bulk material with m layers, respectively. A stands for the in-plane area of the (4) E exf (n) = E iso (n) − E bulk n∕m A bulk unit cell. First, we calculated the exfoliation energy (25.50-28.05 meV/Å 2 ) of n-layer AB-stacked graphene, as presented in Fig. 7. The values are slightly larger than the previous value (21.07-23.62 meV/Å 2 ) [43], but are close to the experimental value (28.75 meV/Å 2 ) [44]. Therefore, our method for calculating exfoliation energy is reliable. And then, we calculated the exfoliation energy of ABstacked FeB 6 bulk from one layer to five layers, as shown in Fig. 7. The obtained values of exfoliation energy are 59.74-70.84 meV/Å 2 , which is significantly smaller than the exfoliation energy value of 2D boron layer separated from bulk KB 4 (150 meV/Å 2 ) [45]. So, the FeB 6 monolayer may be obtained via mechanical stripping. At last, to facilitate comparison with experimental data in the future, the simulated X-ray diffraction (XRD) powder patterns based on AB-stacked FeB 6 bulk were prepared using the Mercury software [46]. First, we calculated the simulated XRD patterns of graphite, as shown in Fig. 8. It can be seen that the graphite has one strong peak (2θ = 26.19°) and two weak peaks (2θ = 42.39°and 44.53°), which consist well with the experimental results [47]. Thus, the method of simulated XRD patterns is reliable. And then, we calculated the simulated XRD patterns of AB-stacked FeB 6 bulk, as shown in Fig. 8. The strongest peak of FeB 6 bulk mainly arises in the range 2θ = 25.82° that is the (1,0,1) reflection plane. And the three weaker peaks are observed at 2θ = 18.89°, 35.31°, and 42.39°, respectively. Besides, the four weakest peaks mainly arise at 2θ = 33.03°, 38.32°, 40.34°, and 49.12°, respectively. These features may be helpful for identifying the FeB 6 bulk in experiments.

Conclusions
A proposal that is the possibility of building graphene analogue FeB 6 with the Fe@B 6 H 6 aggregates is presented theoretically in the present study. We have studied the geometric and electric structures, chemical features, and stability of the monomer, dimer, and trimer of Fe@B 6 H 6 using density functional theory. The results indicate that adopting the dehydrogenation route to polymerization, we can not only place a Fe atom into hexagonal hole instead of B atom in boron sheet, achieving the uniform chemical doping, but also obtain a very large hexagonal-hole density that is very near to 1/3. Moreover, the graphene analogue FeB 6 with 0.4 eV band gap is predicted to be dynamically and thermally stable. The relatively modest exfoliation energy (59.74-70.84 meV/Å 2 ) indicates that the FeB 6 monolayer may be obtained via mechanical stripping. The simulated XRD patterns of FeB 6 bulk provide a reliable prediction for experimental identification. Thus, the new graphene analogue FeB 6 await future experimental verification.
Author contribution Chao Wang: performed the data analyses and wrote the manuscript.
Qiyue Chen: performed the calculations of FeB 6 monolayer. Jianhua Hou: contributed to the conception of the study. Qian Duan: helped perform the analysis with constructive discussions. The n-layer exfoliation energies per area of graphene and FeB 6 . The diagram of the exfoliation process of ABstacked FeB 6 bulk Fig. 8 The simulated XRD patterns of graphite and FeB 6 bulk. The X-ray wavelength is 1.5406 Å