Bonding-restricted structure search for novel 2D materials with dispersed C2 dimers

Currently, the available algorithms for unbiased structure searches are primarily atom-based, where atoms are manipulated as the elementary units, and energy is used as the target function without any restrictions on the bonding of atoms. In fact, in many cases such as nanostructure-assembled materials, the structural units are nanoclusters. We report a study of a bonding-restricted structure search method based on the particle swarm optimization (PSO) for finding the stable structures of two-dimensional (2D) materials containing dispersed C2 dimers rather than individual C atoms. The C2 dimer can be considered as a prototype of nanoclusters. Taking Si-C, B-C and Ti-C systems as test cases, our method combined with density functional theory and phonon calculations uncover new ground state geometrical structures for SiC2, Si2C2, BC2, B2C2, TiC2, and Ti2C2 sheets and their low-lying energy allotropes, as well as their electronic structures. Equally important, this method can be applied to other complex systems even containing f elements and other molecular dimers such as S2, N2, B2 and Si2, where the complex orbital orientations require extensive search for finding the optimal orientations to maximize the bonding with the dimers, predicting new 2D materials beyond MXenes (a family of transition metal carbides or nitrides) and dichalcogenide monolayers.

quite chemically active and inclined to aggregate, it naturally occurs in the carbon vapor of interstellar medium and electric arcs. Vapor from a 16-V carbon arc was found to contain 28 wt% diatomic carbon 13 . As a result, in molecular dynamics (MD) simulations of growth of carbon-based materials, such as graphene, carbon nanotube and C 60 , C 2 dimers are usually adopted as the initial carbon source [14][15][16] .
Furthermore, the C 2 dimer is the building block of metal-alkynide complexes 17 , alkynide complexes 18 , metallocarbohedrenes clusters M m C n 19,20 (usually termed as metcars) and ternary metal carbides (LiAgC 2 , KAgC 2 , CsAgC 2 and NaPdC 2 ) 21 . In recently studied silicon carbides (Si x C y ), the dimerization of C in the SiC monolayer reintroduces the Dirac cone into the honeycomb lattice 22 , while SiC sheets without C 2 dimer are wide-band gap semiconductors [23][24][25][26] . In SiC 2 silagraphene 27 , each silicon atom binds to four C 2 units in a flat plane, resulting in a metallic sheet. When C 2 dimers are embedded in a 2D porphyrin sheet, the system becomes highly active for oxygen reduction reaction 28 . In addition, it was demonstrated that the formation of C 2 dimers plays an important role in stabilizing metcar clusters 29 .
Due to the C 2 dimer's unique properties as well as its importance in forming numerous carbon-based three-dimensional materials and zero-dimensional nanostructures, we focus this study on developing a dimer-based global search algorithm to design 2D materials containing dispersed C 2 dimers to retain its intrinsic properties. We then apply this method to identify the ground state geometries of the 2D Si-C 2 , B-C 2 , and Ti-C 2 systems containing C 2 dimers as test cases.

Results and Discussion
To test our method, we applied it to three systems where C 2 dimers are linked with Si, B and Ti atoms respectively. New ground state geometries and new low-lying isomers are identified.
2D Si-C 2 system. Recently, SiC 2 and Si 2 C 2 sheets with dispersed carbon dimers have received some attention.
The studies were carried out using a priori geometric structures to investigate the properties of 2D Si-C sheets 22,27,30 . However, a complete understanding of Si-C 2 phases is still lacking, and it is unclear whether the studied structures are the ground states under the condition of discrete C 2 dimers. Thus we have conducted a comprehensive structure search to identify the ground state geometry and the low-lying energy isomers of the Si-C 2 system.
The most stable structure of Si 2 C 2 and three low-lying energy isomers of SiC 2 have been identified using our C 2 dimer-based global search. Because the lowest energy geometry of Si 2 C 2 is found to be identical to a recently reported structure 22 , we only concentrate on the three isomers of the SiC 2 sheet, as shown in Fig. 1. In these structures all carbon atoms are pairwisely bonded with each other and covalently bonded with Si atoms. They all are comprised of four-fold coordinated Si with three-fold coordinated C atoms. The structure, shown in Fig. 1(a), is composed of pure pentagons, thus it is named penta-SiC 2 . The other two isomers shown in Fig. 1(b,c) consist of 4-, 5-and 6-rings, and 4-and 6-rings, respectively, labeled as 456-SiC 2 and 46-SiC 2 accordingly. Si-C and C-C bond lengths are about 1.90 Å and 1.36 Å respectively. The side views show the sandwich-like structures with the four-fold coordinated Si atoms in the middle sandwiched between C 2 dimers. When we consider the Si sublattice and C 2 dimer (treating C 2 as a structural unit) sublattice separately, it turned out that the major difference among these three isomers is the orientation of C 2 dimers, namely the dimers are parallel to each other in 46-SiC 2 , while they are perpendicular to each other in penta-SiC 2 (In Supplementary Information, we shortly discuss the kinetic barriers connecting different isomers). It is obvious that the different arrangement of C 2 dimers remarkably alters the geometrical structures, resulting in the different electronic properties of these isomers as demonstrated in the following paragraph.
Total energy calculations reveal that all three isomers are energetically more favorable than the previously proposed SiC 2 -silagraphene 27 with planar tetracoordinate silicon atoms. Penta-SiC 2 , 456-SiC 2 and 46-SiC 2 are 0.63, 0.57, and 0.51 eV/formula unit (f.u.) lower in energy than that of SiC 2 -silagraphene, suggesting that the Si-C 2 sheet is more likely to adopt a buckled structure with the partially sp 3 hybridized Si and the sp 2 hybridized C. Penta-SiC 2 is found to be the lowest energy configuration among the three isomers with the Si-C and C-C bond lengths of 1.90 Å and 1.36 Å respectively, which is identical to that reported recently by Bezanilla et al. 30 , and as 456-SiC 2 lies only 0.06 eV higher in energy than penta-SiC 2 , the two structures can be considered as energetically degenerate. The relative energy of 46-SiC 2 with respect to penta-SiC 2 is 0.12 eV. The dynamical stability of the predicted low-lying energy isomers of SiC 2 is verified by calculating their phonon dispersions. The results plotted in Fig. S1 (see Supplementary Information) show that the three structures are all dynamically stable since their vibration modes are real in their entire Brillouin zones, respectively.
To understand the electronic properties of the three isomers, their band structures are calculated using the HSE06 hybrid functional. Furthermore, a more detailed investigation of the contribution from each atomic orbital to the band structure is carried out using the color-mapped bands. The calculated results are presented in Fig. 1, which shows that penta-SiC 2 , 456-SiC 2 are large band gap compounds with band gaps of 3.1 eV, 2.7 eV, while the bandgap of 46-SiC 2 is 0.55 eV. We note that 46-SiC 2 is a direct band gap semiconductor and its band gap is much smaller than that of penta-SiC 2 and 456-SiC 2 . While the others are indirect band gap semiconductors because the valence band maximum (VBM) and the conduction band minimum (CBM) are not at the same point in their Brillouin Zones. The localized distribution of p z orbitals of the C 2 dimer, above and below the Fermi level, is disclosed in the colored band structure (red dot curves in decomposed band structure). Further analyses of the chemical bonding of the C 2 dimer reveals that the localized distribution originates from the bonding and anti-bonding p z orbital of C 2 dimers. Due to the unique geometrical structure, the alternation of partially sp 3 hybridized Si and slightly distorted sp 2 C prevents the p z orbitals from forming a delocalized π bond. This could explain, to some degree, why these p z composed bands are less dispersive than those in graphene or other planar silicon carbon sheets.
To understand the origin of the different electronic structures exhibited in the three isomers, the detailed structural configurations of the local C and Si atoms in penta-SiC 2 and 46-SiC 2 are sketched in Fig. S3 (see Supplementary Information). Penta-SiC 2 has higher symmetry subjecting it to the least geometrical distortion among the three isomers. 46-SiC 2 's parallel arrangement of the C 2 dimer imposes a larger geometrical distortion; the angle of C with its neighboring Si is 92.9 °, while that in penta-SiC 2 is 109.6 °. The Si atoms in 46-SiC 2 also adopt the distorted hybridization as compared to that in penta-SiC 2 . Besides, the decomposed band structures ( Fig. 1) show that the highest and the second highest occupied bands near the M point in 46-SiC 2 are exclusively from the Si-C σ bonds. The larger geometrical distortion of Si-C σ bonds in 46-SiC 2 leads to a stronger interaction between the relevant orbitals, resulting in the highly dispersive bands and a smaller band gap.
2D B-C 2 system. For the second test system, we study the 2D B-C 2 system due to boron's rich chemistry: for example its electron deficiency and versatile bonding ability. In fact, nanostructures composed of boron and carbon atoms have attracted growing interest both theoretically and experimentally. For instance, a BC 3 honeycomb structure with high crystalline quality was identified by electron diffraction 31 . A B 2 C sheet with planar tetra-coordinated carbon (pt-C) moiety was theoretically predicted 32 . A thorough study of B-C compounds with different stoichiometry was carried out by performing a PSO search 33 . However, a study on B-C sheets comprised of discrete C 2 dimers has not been reported. Due to its versatility in forming various chemical bonds and the differing chemical nature of B from Si, the bonding between B and the C 2 dimer is expected to be quite different from that of Si and C 2 , leading to B-C 2 systems with different structural and electronic properties.
Using our dimer-based search algorithm, three low-lying energy isomers of B 2 C 2 and one stable structure of BC 2 are obtained (Fig. 2). Geometry optimizations and total energy calculations suggest that the structure shown in Fig. 2(a1) has the lowest energy of the three isomers of B 2 C 2 with a total energy of 0.04 eV, 0.13 eV lower than that of the structures shown in Fig. 2(a2,a3), respectively. To examine the dynamical stability of these structures, we calculated their phonon dispersions. The results are plotted in Fig. S4. The absence of imaginary modes in the entire Brillouin zone for each of the structures confirms that they are all dynamically stable. However, we note that these structures are energetically metastable as compared to the B-C planar monolayers predicted by Luo et al. (B 2 C 2 -a is 0.23 eV/f.u. higher in energy than the most stable BC sheet).
In order to gain more insight into the chemical bonding in the B-C 2 sheets, we plotted an electron localization function (ELF) 34 isosurface (with a relatively large isovalue of 0.75 to reflect the σ bonds between the atoms) in Fig. 3. According to the ELF diagrams, carbon in B 2 C 2 -a and B 2 C 2 -b adopt sp 3 hybridization, and each carbon atom is covalently bonded to three adjacent B and one C atoms. B atoms take a slightly distorted sp 2 configuration. The C-C bond length in Fig. 2(a1,2) are 1.51 Å (close to that of a C-C single bond). While in B 2 C 2 -c, both the B and C atoms have sp 2 hybridization. B 2 C 2 -c adopts a planar honeycomb structure similar to graphene 35 , however, the total energy of B 2 C 2 -c is 0.13 eV/f.u. higher than that of B 2 C 2 -a. Multi-center bonds are pretty common in boron-related materials 36,37 , and contribute a lot to their stability. However, no such bonds are found in B 2 C 2 -a, B 2 C 2 -b and B 2 C 2 -c, which might explain the physical origin of their higher energy.  For the BC 2 sheet, only one stable structure without imaginary modes was discovered (Figs 2(b) and S4(b)). In particular, all the BC 2 structures isostructural to those of SiC 2 (namely penta-BC 2 or 46-BC 2 ) are unstable, reflecting the differing chemical nature of B and Si. In this highly buckled BC 2 sheet (the thickness is 1.53 Å), all C atoms are strongly covalently bonded together; the C-C bond length is 1.31 Å. The calculated electron localization function (ELF) shows that there are sp 3 hybridized B and nearly planar tetracoordinate B atoms. In our case, the bonding type of planar tetracoordinated B is closer to the case of planar tetracoordinated Si 27 and C 32 .
Our electronic band structure calculations (Fig. 2) show that the ground state structures of B 2 C 2 and BC 2 with the discrete C 2 dimers have indirect band gaps of 2.9, 2.8 and 4.1 eV at the HSE06 level for B 2 C 2 -a, B 2 C 2 -b and BC 2 , respectively, which is in contrast to the metallic BC sheet 33 where C is not in C 2 form. The reason is, once we restrict the carbon atoms to be pairwisely bonded together, the most stable structures (B 2 C 2 -a and BC 2 ) discovered by using our dimer-based search are buckled. The separation of C 2 dimers by boron atoms and the buckling nature of the sheets would destroy the delocalization of p z π electrons of boron and carbon (the main reason for the metallicity of many planar stable B-C compounds). In the B 2 C 2 -c isomer, the planar geometry, as well as its delocalized π electrons is preserved. Geometrically speaking, B 2 C 2 -c is similar to honeycomb SiC 22 . Since boron has few valence electrons it is impossible to fill all the electronic states below the Dirac cone giving rise to metallicity.
2D Ti-C 2 system. To study the interaction between the C 2 dimer and transition metal elements, we have systematically studied Ti-C 2 sheets for the following reasons: (1) both Ti and Si possess four valence electrons, thus it is intriguing to look at the difference between TiC 2 and SiC 2 sheets and Ti 2 C 2 and Si 2 C 2 sheets. (2) Titanium and carbon have been found to form various kinds of stable compounds ranging from 0D clusters 19,38 , to 2D sheets 39,40 , and to 3D bulk materials. Among them, in metcar clusters all carbon atoms exist in the form of a C 2 dimer 19 . Therefore, it might be possible to find a family of 2D crystals that consists of transition metal atoms and C 2 dimers. (3) Recently, a group of 2D layered early transition metal carbides termed as MXenes, including Ti 3 C 2 and Ti 2 C, were experimentally synthesized by exfoliation from their bulk MAX phases [41][42][43] . However, in MXenes sheets, C is in atomic form rather than a C 2 dimer. Therefore it is necessary to study the possibility of forming stable 2D titanium carbides containing C atoms that are all in C 2 dimer form.
We concentrated on the 2D titanium carbide sheets with the stoichiometry of TiC 2 and Ti 2 C 2 where C is in the form of a C 2 dimer. Recently, a TiC 2 sheet was predicted to have good performance as anode material for a lithium ion battery 40 . However, the geometrical structure was artificially designed. Thus, it is easy to miss some low-lying isomers. In fact, we found, using our dimer-based algorithm, that among the isomers of TiC 2 , there is a one that is lower in energy by 0.13 eV/f.u. than the recently reported structure 40 . In this lowest energy structure, as shown in Fig. 4(a), due to the different orientations of C 2 dimers, the total thickness of 1.97 Å is less than the 2.02 Å of the reported structure 40 . One can see when going from Si to Ti, although both have four valence electrons, the complex 3d orbitals in Ti require a delicate arrangement of C 2 dimers to match with the d orbital orientations for maximizing the interactions. The stability of the TiC 2 sheet can be traced back to the most stable C v 3 -like Ti 8 C 12 cluster 44 , where all carbon dimers bind to Ti in end-on configuration (EOC) and side-on configuration (SOC), which strengthens the bonding between Ti and C and stabilizes the Ti 8 C 12 cluster. In our 2D TiC 2 sheet EOC and SOC are apparently well preserved. The C-C bond length is 1.33 Å, which is slightly longer than that of a dispersed C 2 molecule (1.31 Å) due to the partial occupation of π g orbital of C 2 in TiC 2 . The bond length of Ti-C in EOC and SOC modes is 2.05 Å and 2.19 Å, respectively, close to 2.1 Å in Ti 2 C MXene phase.
We then increased the atomic ratio of Ti from TiC 2 to Ti 2 C 2 to study the effect of Ti concentration on the geometry and properties of the Ti-C sheets. The structure shown in Fig. 4(b) is found to have the lowest energy among the isomers of Ti 2 C 2 under the requirement of discrete C 2 dimers, and it lies lower in energy by 0.25 eV/f.u than the t-TiC sheet 39 . This suggests that when the same atomic ratio of Ti to C exists, C 2 is energetically more favorable than atomic C to bind with Ti in forming a sheet. The C-C bond length is 1.48 Å, close to that of C-C single bond but longer than that in a TiC 2 sheet. This is because the hybridization between the C 2 -π g and Ti-3d is stronger in Ti 2 C 2 , resulting in more π g states being occupied, thus leading to the elongated C-C bond.
The calculated band structures and partial DOS are also plotted in Fig. 4, which clearly show that the TiC 2 sheet is metallic, while the Ti 2 C 2 sheet, like the t-TiC sheet, is an indirect narrow gap semiconductor with a gap of 0.1 eV at the HSE06 level. It is known that the ground state electronic configuration of an dispersed C 2 is (1σ g ) 2 (1σ u ) 2 (2σ g ) 2 (2σ u ) 2 (1π u ) 4 with the higher energy states (3σ g )(1π g )(3σ u ) unoccupied, and the (2σ g ) 2 (2σ u ) 2 (1π u ) 4 (3σ g ) having much lower energy as compared to the orbitals of Ti-3d. Therefore, a weak interaction between the deep level orbitals of C 2 and Ti-3d is expected. Thus, the electronic states near the Fermi level are mainly from the interplay between C 2 -π g and Ti-3d orbitals.
For the structure of t-TiC 39 a strong electronic resonance between C-2p and Ti-3d is obvious. Due to the strong hybridization between C-2p and Ti-3d and the low buckled geometry of t-TiC, the bonding states of C-2p and Ti-3d are separated from the unoccupied Ti-3d states forming a band gap. TiC 2 is geometrically analogous to t-TiC with carbon atoms replaced with discrete C 2 dimers. However, replacing C with C 2 leads to drastic change in their electronic structures, namely t-TiC is semiconducting with a strong electronic resonance, while TiC 2 is metallic. This is because the π g orbitals of C 2 hybridize with the Ti-3d orbitals, forming the energy bands crossing the Fermi level, thus, resulting in metallicity of TiC 2 . For the Ti 2 C 2 sheet, the electronic states near the Fermi level are also derived from C 2 π g and Ti 3d orbitals similar to TiC 2 , yet the interaction between C 2 π g and Ti 3d orbitals is stronger as compared to that in TiC 2 due to the higher Ti concentration, which reduces the electronic states near the Fermi level, and thus results in the small band gap as shown in Fig. 4(b).
Based on the above discussions, the following intriguing points can be concluded: (1) compared with the sandwiched structure of the TiC 2 sheet, the Ti atoms are more exposed in the Ti 2 C 2 sheet, and hence have higher chemical reactivity; (2) Although having a higher atomic ratio of Ti, the Ti 2 C 2 sheet is semiconducting with a band gap of 0.1 eV, while the TiC 2 sheet is metallic; (3) For the t-TiC 39 and Ti 2 C 2 sheets, although they have the Scientific RepoRts | 6:29531 | DOI: 10.1038/srep29531 same composition ratio, when going from atomic C to C 2 dimer, the band gap is reduced from 0.2 eV to 0.1 eV, showing promising applications of the Ti 2 C 2 sheet in infrared light-related technology.

Summary.
In this study, we developed a dimer-based bonding-restricted structure search method to find the ground state structures of the 2D materials containing C 2 dimers with a restriction of no direct bonded dimers, and then applied this method to three systems composed of discrete C 2 dimers that are coordinated with Si, B and Ti atoms, respectively, as test cases. For the SiC 2 sheet, three energetically nearly degenerate allotropes with very different electronic structures are identified, which also differ from the previously proposed SiC 2 -silagraphene 27 , showing that tuning the orientation of the C 2 dimers not only can modulate the size of the band gap, but also can induce the transition from an indirect gap to a direct gap. For the B-C system, the ground state structures of both BC 2 and B 2 C 2 sheets are determined. Due to the different chemical nature of B and Si, all the possible structures of BC 2 isostructural to those of SiC 2 are unstable. In the ground state configuration of BC 2 , each C 2 dimer is bonded with three B atoms forming a 2D sheet with a thickness is 1.53 Å and a band gap of 4.1 eV. When increasing the number of atomic B to form B 2 C 2 , each C 2 dimer is bonded to six B atoms having the maximum number of B-C σ bonds and reducing the band gap to 2.9 eV. Unlike the metallic BC sheet 33 where C atoms do not form any dimers, both the BC 2 and B 2 C 2 sheets identified here are semiconducting. The emergence of band gap is attributed to the buckled structures, which prevents the delocalization of p z orbitals. For the Ti-C system, a new structure of the TiC 2 sheet with a lower energy than the metallic TiC 2 sheet reported recently 40 was found. A Ti 2 C 2 sheet was also found to be energetically more stable than the previously proposed t-TiC sheet 39 , suggesting it is more favorable for Ti to bind with C 2 dimers in forming 2D structures. Unlike the metallic TiC 2 monolayer 40 , the Ti 2 C 2 sheet is an indirect band gap semiconductor. In addition, compared with the structures composed of C 2 dimers and nonmetallic elements (Si, B), both the stable TiC 2 and Ti 2 C 2 sheets adopt relatively more compact structures due to the complex orientations of d orbitals, which require a more extensive search to find the optimal orientation of the C 2 dimers. These cases show that our searching algorithm is efficient and indispensable to design new 2D materials beyond the atom-based ones.

Methods
Dimer-based search algorithm. Particle Swarm Optimization (PSO), one of the most popular swarm intelligence algorithms, was originally developed by Kennedy and Eberhart in 1995 45 . The PSO algorithm utilizes the collective intelligence of the whole population rather than a single individual, it is quite efficient in handling various optimization problems. Recently, the PSO algorithm has been successfully applied to crystal structure prediction 46-48 , another intractable optimization problem. Just like in the natural world where collective intelligence facilitates the locating of food or optimal habitats for insects and fishes, PSO implementation in CALYPSO was proven to be efficient in predicting the most stable structures for undiscovered materials 48 . Due to the great success of PSO in predicting crystal structures, we have studied a dimer-based structure search program based on the PSO algorithm, where additional constraints and augmented algorithm are added. The main framework of our dimer-based structure search method is illustrated by the flowchart as shown in Fig. 5.
Before outlining the details of the steps given in the flowchart, we will explain the treatment of the C 2 dimer in our constrained search method. The requirements are that all carbon atoms exist in the form of a C 2 dimer and each dimer must be separated from the others. A C 2 dimer can be treated naturally as a "pseudo-atom". However, general treatments for ordinary atoms are not sufficient for "pseudo atoms", as an ordinary atom is treated as an isotropic sphere, while a C 2 dimer is considered a pseudo atom exhibiting structural anisotropy along the molecular axis. Therefore, to fully represent a C 2 dimer, extra indices are imperative. As shown in Fig. 6, besides the three coordinates, x, y and z, the dimeric orientation angles ϕ and θ are required. Therefore, in our augmented PSO algorithm, the coordinates needed are r i = (x i , y i , z i , ϕ i , θ i ), where x i , y i and z i are used to denote the barycenter, ϕ i and θ i are used to represent the orientation of the dimer. Besides, during the PSO-based structure searches, the bond length of C 2 dimer is allowed to vary within the length of C-C single and triple bond. Accordingly, the PSO algorithm is modified to accommodate these changes.
There are four main steps enumerated in the flowchart, some operative techniques are directly adapted from previous works 46,48 . First of all, most of the trial structures are partially randomly generated under symmetry restrictions 46,49,50 or evolved via the PSO algorithm 48,51 . Fully randomly generated structures are also used. The reason is in 2D space there are 17 plane groups overall belonging to four different crystallographic systems (oblique, rectangular, square, and hexagonal); one of the groups is selected each time to create a new trial structure. The lattice shape of the new structure is determined by the selected space group and its area. Generally, the positions of the atoms or the barycenter of C 2 dimer are allocated on Wyckoff positions for the chosen plane group. An extra step is needed to generate bond length and orientation for each dimer randomly or partially restricted, which is different from atom-based structure prediction methods. All these steps are sketched in Fig. 7. In order to predict 2D materials with finite thickness, the coordinate z i of atoms and dimers are assigned randomly within a thickness Δ as shown in Fig. 7(d). However, because the number of 2D plane groups is finite, the generating method for a symmetry constrained structure is prone to produce similar structures. In this case totally randomly generated structures are needed to guarantee the diversity of structures.
Another source of the new trial structures is the PSO evolution which is a population-based algorithm. The population in our case is defined as a group of structures. For generation t, we define a population P t and every structure in generation t is represented as a particle P i t . The evolution of structures is controlled by PSO algorithm: where i and t are indexes of particle (structure) and generation, R i t represents the structural information of a certain structure including the orientation of dimers and the coordinates of atoms, v i t is the velocity or rate of structural change of certain structure i. The velocity of structure i is calculated as follows: where ω is the inertia weight and controls the momentum of the particle, R pbest t i is the optimized structure of particle i, and R gbest t i is the structure with the lowest energy obtained so far for generation t, r 1 and r 2 are two random numbers uniformly distributed in range [0, 1], c 1 and c 2 are set to be 2 47,52 . The definition of velocity fully reflects the idea of collective intelligence. The trajectory of each particle is biased, to a certain degree, toward the best structure of whole population and itself at the same time.  The three structure generation approaches work synergistically and serve different purposes during the structure search. The partially random structure generation with symmetry constraints help reduce the degree of freedom and narrow the search space. The fully random structure generation warrants the structural diversity, and the PSO-based structure evolution provides an effective way to explore the configuration space and locate the thermodynamically favorable structures.
It is important to note here that due to the nonequivalent coordinates x i , y i , z i , and ϕ i , θ i in terms of their magnitudes and range, C 2 dimers and ordinary atoms must be manipulated separately. In addition, the newly generated trial structures must satisfy the condition that all C 2 dimers are disconnected from each other, thus the structures with dimer clustering are eliminated. Next, similarities between the different structures are examined. Each structure is characterized by a set of predetermined fingerprinting functions 48,52 , and the distance between two structures are calculated based on the difference of their fingerprinting functions. If a newly produced structure is similar to the previous one, the new structure is discarded. Finally, local structure optimization is performed to drive the structure to the local minimum on multi-dimensional energy space, which also provides physical and rational structures for future PSO evolution.
First principles calculation. Combined with our bonding-restricted algorithm, first principles calculations within the framework of density functional theory (DFT) are performed for local geometry optimizations during structure search. Calculations are mainly performed by using the Vienna Ab initio Simulation Package (VASP) 53 . The Perdew-Burke-Ernzerhof (PBE) functional 54 is used to incorporate the exchange-correlation interaction. The HSE06 hybrid functional 55,56 is used for high accuracy of electronic structure calculations. The projector augmented wave (PAW) 57 method is used to treat the interactions between ion cores and valance electrons. Plane waves are used to expand the valance electron (2s 2 2p 1 for B, 2s 2 2p 2 for C, 3s 2 3p 2 for Si, and 3d 3 4s 1 for Ti) wave functions. During the massive structure searching stage, in order to reduce the workload, plane waves with a kinetic energy cutoff of 350 eV and Monkhorst-Pack scheme with a sparse grid density (2π × 0.04 Å −1 ) are adopted. To compare the relative stability of different candidate structures, we used the same kinetic energy cutoff, k-point grid density and exchange-correlation functionals to perform geometry optimizations and total energy calculations. Plane waves with a kinetic energy cutoff of 450 eV and Monkhorst-Pack scheme 58 with a grid density of 2π × 0.02 Å −1 are used to optimize the structures and calculate their electronic properties. The convergence criteria for total energy and forces are set to be 10 −4 eV and 10 −3 eV/Å, respectively. The supercell approach is used to calculate force constants. A vacuum space of 20 Å in the perpendicular direction is added to avoid the interaction between periodic images. Phonon properties are calculated using the finite displacement method as implemented in the Phonopy package 59,60 .