Prediction of unexpected BnPn structures: promising materials for non-linear optical devices and photocatalytic activities

In the present work, a modern method of crystal structure prediction, namely USPEX conjugated with density functional theory (DFT) calculations, was used to predict the new stable structures of BnPn (n = 12, 24) clusters. Since B12N12 and B24N24 fullerenes have been synthesized experimentally, it motivated us to explore the structural prediction of B12P12 and B24P24 clusters. All new structures were predicted to be energetically favorable with negative binding energy in the range from −4.7 to −4.8 eV per atom, suggesting good experimental feasibility for the synthesis of these structures. Our search for the most stable structure of BnPn clusters led us to classify the predicted structures into two completely distinct structures such as α-BnPn and β-BnPn phases. In α-BnPn, each phosphorus atom is doped into a boron atom, whereas B atoms form a Bn unit. On the other hand, each boron atom in the β-phase was bonded to a phosphorus atom to make a fullerene-like cage structure. Besides, theoretical simulations determined that α-BnPn structures, especially α-B24P24, show superior oxidation resistance and also, both α-BnPn and β-BnPn exhibit better thermal stability; the upper limit temperature that structures can tolerance is 900 K. The electronic properties of new compounds illustrate a higher degree of absorption in the UV and visible-region with the absorption coefficient larger than 105 cm−1, which suggests a wide range of opportunities for advanced optoelectronic applications. The β-BnPn phase has suitable band alignments in the visible-light excitation region, which will produce enhanced photocatalytic activities. On the other hand, α-BnPn structures with modest band gap exhibit large second hyperpolarizability, which are anticipated to have excellent potential as second-order non-linear optical (NLO) materials.


Introduction
In the past decade, many attempts have been made to discover non-carbon fullerene-like structures based on their unique properties. [1][2][3][4][5][6][7][8][9] The structural and electronic properties of clusters form an intermediate state between the individual atoms and crystals. Clusters have a remarkable impression in understanding the chemical bonding and the rational design of molecules with tailored physicochemical properties. In recent years, scientists have focused on studying boron-based materials because a signicant number of unfamiliar structures of boron and its compounds have been explored. [10][11][12][13][14][15][16] The electron conguration of boron is 2s 2 2p 1 ; therefore, boron with one electron less than carbon has different characteristics compared to carbon, which can form fascinating and strong bonds with nearly all periodic table elements. [16][17][18] On the other hand, phosphorus clusters with large pore sizes have received much attention in the past decade. 18 Phosphorus is the diagonal neighbor of carbon in the periodic table; therefore, the chemical bonding and structural chemistry of phosphorus are also interesting. Elemental phosphorus exists in different forms such as white, yellow, and black phosphorus but phosphorus is never found as a free element on earth because of its high reactivity. 19 The group III-V compounds such as boron phosphide (BP) have been most generally investigated in their cubic crystal phase (c-BP). 20 Also, BP compounds can exist in planar (twodimensional), tubular (one-dimensional), or spherical shapes (zero-dimensional), quite the same as graphene, carbon nanotubes, and fullerenes, respectively. All structures of BP compounds are prominent as large band-gap semiconductor materials having many advanced applications including microelectronics and optoelectronics, 21,22 high thermal conductivity 23 hardness, chemical stability, [23][24][25] and catalysis. 20,26,27 Both boron and phosphorus elements can form 2D materials; it is interesting to determine whether the binary compound of boron and phosphorus can form stable clusters, which may display unusual structural and electronic properties better than both boron and phosphorus clusters. 28 Hence, it is desirable to consider the structures of boron phosphide clusters with different sizes.
Numerous computational and experimental studies have considered the structural and electronic properties of group III-V compounds. 18,19,[29][30][31][32] In recent years, non-carbon (inorganic) (XY) n fullerenes of group III-V are of great scientic interest based on their promising candidates as optoelectronic devices 33,34 and light-emitting diodes (LEDs). 35 X 12 Y 12 , X 16 Y 16 , X 24 Y 24 , and X 28 Y 28 clusters, where X ¼ group III atom and Y ¼ group V atom, are predominant inorganic (XY) n fullerenes. 8 Computational investigations on (XY) n fullerenes have demonstrated that (XY) 12 has the highest thermodynamic stability, which is known as a magic cluster. [36][37][38] In this regard, Strout et al. ascertained that the X 12 Y 12 nano-cages are the most stable structures between all fullerene-like (XY) n (X ¼ B, Al, Ga,. and Y ¼ N, P, As,.) structures. 36 B 12 N 12 and B 24 N 24 fullerene-like clusters were synthesized by Oku et al. [39][40][41] and detected using time-of-ight mass spectrometry. The B 12 N 12 cluster is shown as the most stable structure among B 12 N 12 , B 16 N 16 , and B 28 N 28 . 42 Other (XY) 12 clusters such as Al 12 N 12 , Al 12 P 12 , B 12 N 12 , and B 12 P 12 are of much interest owing to their unique chemical and physical properties. In this regard, many theoretical studies had been focused on ascertaining the relative stabilities of different sizes of these nano-cages. 18,19,[30][31][32]38,43 Boron phosphide (BP) n clusters are an attractive member of the X n Y n nano-cage family with unique physical and electrochemical properties. Since B 12 N 12 and B 24 N 24 fullerenes have been synthesized experimentally, [39][40][41] it has motivated researchers to investigate the structural and electronic properties of B 12 P 12 and B 24 P 24 fullerene-like clusters. In all the previous studies, X 12 Y 12 and X 24 Y 24 clusters have been considered as fullerene-like structures, which consist of four-, six-, and eight-membered rings, and have large HOMO-LUMO gaps, almost zero dipole moment, and nearly zero hyperpolarizability. 31 It should be mentioned that in all the previous studies, (XY) n clusters (X ¼ B, Al, Ga, and Y ¼ P, As) were made directly by replacing B and N atoms of the already known structures of the (BN) n cages with other atoms, then the bonding length and angle are adjusted by further calculations. Based on this strategy, the main questions are that whether this method is correct and whether this method give the lowest energy structure. We believe that structural prediction can answer these questions. Structural prediction based on rstprinciples calculation has served as a very useful tool in materials research. Exploring 0, 1, and 2-dimensional (0, 1, 2-D) materials from molecular design and global search have been a hot-topic. Based on our knowledge, no study has been conducted to predict the structure of boron phosphide (BP) n clusters. Also, all previous studies on these materials have been based on the structure of boron nitride (BN) n clusters. The structure of B n P n in all previous investigations by other researchers has been reported to be cage-shaped, which is exactly similar to the (BN) n clusters, in which the P atoms are replaced by N atoms. It motivated us to use an evolutionary algorithm to predict the structures of boron phosphide clusters. For the rst time, we used the evolutionary algorithm as implemented in the USPEX code to search the ground state structure of the B n P n clusters, including B 12 P 12 and B 24 P 24 clusters. It should be noted that the B 12 N 12 and B 24 N 24 fullerenes have been synthesized experimentally; hence, the B 12 P 12 and B 24 P 24 structures were chosen.

Computational details
We used the evolutionary algorithm as executed in the USPEX code to search for the potential energy surface of the B 12 P 12 and B 24 P 24 clusters. USPEX code by evolutionary algorithm solves the central problem of crystal structure prediction of theoretical crystal chemistry. [44][45][46][47] An adequate number of initial structures, more than 1800 congurations in 50 generations by the USPEX method, was used to ensure that the global minimum is attained. The geometrical relaxation, stability assessment, and optoelectronic properties of the predicted B n P n clusters were carried out by the Vienna ab initio simulation package (VASP). 48,49 The generalized gradient approximation (GGA) scheme, developed by Perdew-Burke-Ernzerhof (PBE) and projectoraugmented wave (PAW) 50 potentials were used. The dispersioncorrected DFT-D3 method with Beack-Jansoon damping was used in all the structure optimizations to account for the weak van der Waals interactions. 51,52 The low-lying energy isomers, whose relative energies are in the range of 0-2.5 eV with respect to the most stable structure, were carefully re-optimized to distinguish the correct global minimum structure of each B n P n cluster. The screen hybrid HSE06 method 53 with the mixing parameter alpha value of 0.25 was carried out to consider the electronic properties and optical absorption of the predicted clusters. An energy cutoff of 400 eV was used for the plane-wave expansion. All geometries were relaxed until the convergence criteria of energy (10 À6 eV) and force on each atom (0.02 eV A À1 ) were satised. It should be pointed out that the structures were located in the middle of a box with a vacuum spacing of $10 A in all the directions to neglect the interactions between the adjacent cells. To further examine the thermal stability proposed for the B n P n clusters, ab initio molecular dynamics (AIMD) simulations using the canonical (NVT) ensemble were performed. The simulations are accomplished using a Nosé-Hoover thermostat at nite 300, 600, and 900 K temperatures for 5000 fs with a time step of 1.0 fs.

Geometry and stability
Based on the USPEX method conjugated with spin-polarized rst-principles calculations, a comprehensive study to explore their ground state structures of B n P n (n ¼ 12, 24) and low-lying isomers of these clusters was performed employing the evolutionary algorithm. The USPEX method has been successful used before to explore a wide range of materials. 54-56 Among the structure search by evolutionary algorithm, more than 1800 different structures for each cluster were predicted. We mainly focused on the low-lying isomers in the energy range of 0-2.5 eV with respect to the most stable structure. We reoptimized these structures to get the precise total energy and the atomic structures. The lowest-energy structures and lowlying energy isomers of the B n P n (n ¼ 12, 24) clusters are displayed in Fig. 1 and 2.
Our search for the most stable structure of B n P n clusters led us to classify the predicted structures into two completely distinct structures. Interestingly, it can be classied into two aand b-phases. In the a-phases, each phosphorus atom is doped into a boron atom while the B atoms form a B n unit. This type of predicted structure named a-B n P n structure, is reported for the rst time. On the other hand, in the b-phase, each boron atom is bonded to each phosphorus atom to make a cage cluster, which is named the b-B n P n structure. The nal relaxed geometry of the a-B n P n and b-B n P n systems with their symmetries, their lowlying energy isomers, and relative stabilities are displayed in Fig. 1 and 2. To assist comprehension, two different views of the ground-state structures, as well as the rst low-lying energy isomers of a-B n P n and b-B n P n are shown in Fig. S1 and S2, † respectively.
Based on our global structure search, the most stable structure of the B 12 P 12 cluster with D 2 symmetry is alpha rhombohedral-shaped, named as a-B 12 P 12 . In this structure, the arrangement of 12 boron atoms in the alpha-tetragonal boron structure made a B 12 building block, while each phosphorus atom bonded to each boron atom of the alpha tetragonal B 12 unit. To date, this is the rst time that this type of structure has been reported for the B 12 P 12 structure. In this beautifully predicted structure, each P atom is covalently bonded to each boron atom of the B 12 unit with an average B-P bond length of $1.98 A. In other words, the surface of the B 12 unit is covered with phosphorus atoms. The average B-B bond lengths in the B 12 unit of a-B 12 P 12 ($1.77 A) is nearly the same as that of the average B-B bond in the B 12 icosahedron (1.75 A). 57 As seen in Fig. 1, each phosphorus atom is covalently bonded to two neighboring phosphorus atoms to make a three-membered ring, in which the phosphorus atoms are lying in a plane. It can be imagined that the triangular rings of phosphorus have grown on the B 12 unit, while each of these P atoms are bonded to each of the B atoms. In this situation, the hybridization of each P atom is sp 3 The calculated average P-P bond length and the P-P-P bond angles of a-B 12 P 12 are about 2.36 A and 60 , respectively, which is comparable to the P-P bond length (2.23 A) in phosphorene. 58 The experimental evidence for our predicted a-B 12 P 12 structure is the alpha tetragonal allotrope of boron determined in 1951 by Hoard et al. 59 The same as a-B 12 P 12 structure, this structure is made up of just one building block of the alpha-tetragonal B 12 boron structure. Another evidence for the new predicted structure is B 13 P 2 . 60 The crystal structure of B 13 P 2 shows that this is an alpha rhombohedral boron structure with two P atoms in the unit cell. Interestingly, B 13 P 2 was made up of just one B 13 building block, in which P atoms connected each building block. 60 The rst low-lying energy structure of B 12 P 12 is fullerene-like with a large cavity. This type of structure was named as the b-B 12 P 12 structure. The structure of b-B 12 P 12 is a polyhedral boronphosphide structure composed of four tetragons and hexagons. As seen in Fig. 1, the stability of the low-lying energy b-B 12 P 12 structure is less than the a-B 12 P 12 conguration with a large amount of energy of about 1.013 eV. Based on our calculations, it has an average diameter of 0.586 nm, while the C 60 fullerene has an average diameter of $0.7 nm. 61,62 Due to the large cavity of b-B 12 P 12 , the same as that of C 60 fullerene, it can be suitable for various applications such as gas storage, encapsulation of atoms, or molecules inside its cage, and endohedral metallofullerenes. As we know, the endohedral metallofullerenes have attracted much attention due to their unique potential applications such as superconductors and non-linear optical (NLO) devices. The relaxed geometry of the b-B 12 P 12 cluster indicated that only B-P bonds with an average bond length of $1.911 A were found in this nano-cage. It should be noted that two types of B-P bonds were found in this structure; the rst is the B-P bond in the hexagonal ring, which is named I (B-P), and the second is the B-P bond in the tetragonal rings, which is named II (B-P). These results are collected in Table 1 and are shown in Fig. S1. †  Table 1 The PBE-D3 structural parameters of the relaxed geometries of the most stable structures and the first low-lying energy isomers of the (PB) n clusters (atom numbering is in accordance with see ESI Fig  The nal relaxed geometry of the ground-state structure and the low-lying isomers of the B 24 P 24 cluster, together with their symmetries and relative stabilities, are shown in Fig. 2 and S2. † By looking at the predicted structures of B 24 P 24 , it was found that the ground state structure of B 24 P 24 is a-B 24 P 24 with a high symmetry of D 3d , in which all the boron atoms are arranged as B 24 unit, while the phosphorus atoms cover the B 24 surface. It should be pointed out that the B 24 motif in a-B 24 P 24 is composed of two hexagons on the bottom and the top, six hexagons surrounding the waist, and 12 trigons alternatively connected to each hexagon. To date, this has been the rst time that this type of structure has been reported for the P 24 B 24 cluster. In this unique predicted structure, each P atom is bonded to one or more boron atoms of the B 24 motif as the surface of the B 24 unit is covered with phosphorus atoms. The average B-P bond length is obtained as $1.99 A, which is comparable with the B-P bonds in phosphinoboranes and related compounds. 63 In detail, two types of B-P bonds such as types I and II can be found in the a-B 24 P 24 structure (Fig. S2 †), and the corresponding values of these bond lengths are collected in Table 1. On the other hand, the calculated values of the P-P and B-B bond lengths are summarized in Table 1. The average P-P bond length is obtained as 2.231 A while for the B-B bonds, we found two different types of B-B in sh3 a-B 24 P 24 structure. I (B-B) h represents the bonded B atoms to each other in hexagons (1.746 A) and II (B-B) T shows the bond length between the boron atoms in trigonal geometry with a bond length of about 1.788 A. These results are comparable with the average-B-B bond in the B 12 icosahedra (1.75 A). 57 For simplicity of visualization, two different views of the B 24 unit without P atoms are shown in Fig. S3. † This result is contradictory to that reported for pristine B 24 clusters, where the quasi-planar and rather irregular polyhedrons are prevalent. 64 It should be reminded that the ionic quasi-planar B 24 structure has been experimentally synthesized. 64 The evidence for our predicted a-B 24 P 24 cluster is the encapsulation of a transition metal in the fullerene-like boron cages. [65][66][67] Yanming Ma and coworkers 65 predicted transition metal-doped B 24 clusters using rstprinciples swarm-intelligence-based structure searches. They found that the low-lying energy structures were generally a cagelike structure. It can be concluded that doping is a signicant factor in determining the geometry of boron-based materials. Fig. 2 represents the selected low-lying energy B 24 P 24 structures with their symmetry and energy relative to the groundstate D 3d a-B 24 P 24 structure. The rst low-lying energy structure belongs to the b-B 24 P 24 fullerene-like structure. As seen from this gure, the b-B 24 P 24 structure with D 4 symmetry is less stable than the ground-state a-B 24 P 24 structure by a large amount of energy ($0.68 eV). It should be reminded that the fullerene-like structure is proposed as the ground-state structure of the X 24 Y 24 clusters. The investigated structures of X 24 Y 24 in all the previous studies 19, [30][31][32][33][34]38,43 are based on the strategy in which B and N atoms are directly replaced with X and Y atoms in the already known structures of the (BN) n nano-cages, and then the bonding length and angle are adjusted for further calculation. In the case of B n P n , our global search indicated that the B n P n fullerene-like structure is a low-lying isomer. The b-B 24 P 24 fullerene-like structure with a 0.859 nm cavity, which is larger than the C 60 cavity (0.7 nm), includes 12 tetragonal, eight hexagonal, and six octagonal BP rings. Table 1 summarizes different boron-phosphorus (B-P) bond lengths of the most stable structures and the rst low-lying energy structures of B 24 P 24 . As can be seen in Fig. 2 and Table 1, two kinds of B-P bonds are recognized for the b-B 24 P 24 structures with B-P bonds between tetragonal/hexagonal rings, named type I (B-P), with an average bond length of about 1.918 A, B-P bonds between tetragonal/octagonal rings, named as type II (B-P) ($1.905 A), and the B-P bonds between hexagonal/octagonal ring, named type III, with bond length of about 1.871 A.
The energy difference between the most stable structure and other competing structural isomers of each compound are essential to gure out the possibility of their synthesis under different experimental conditions. The energy difference between the lowest energy structure of a-B 12 P 12 /a-B 24 P 24 and the rst low-lying energy b-B 12 P 12 /b-B 24 P 24 structures is 0.059/ 0.061 eV per atom, respectively, implying that the a-B n P n structures are more stable than the b-B n P n structures ($5.7 kJ mol À1 ). Based on our predictions, b-B 12 P 12 and b-B 24 P 24 clusters may not coexist with the a-B n P n structures at room temperature when compared with the thermal energy at 300 K (0.026 eV). Hence, it can be concluded that the B n P n fullerene-like structures can exist at higher temperatures (>300 K). In the case of other low-lying energy isomers of B 12 P 12 and B 24 P 24 structures, the second and third metastable structures are less stable than the global minimum a-B 12 P 24 and a-B 24 P 24 structures by a large amount of energy (>2.2 eV). Now, we return to investigate the thermodynamic stability of the a-B n P n and b-B n P n clusters, which need to certify their stability and facility experimental synthesis. To estimate the thermodynamic stability of the predicted clusters, the binding energy per atom (E b , eV per atom) is calculated. The relative binding energy (E b ) is an effective parameter to show the thermodynamic stability of the clusters and it is dened as below.
where E (BP) n and E B /E P illustrate the total energies of the B n P n clusters and a single boron/phosphorus atom, respectively. As seen in Table 2, the calculated binding energies of the a-B n P n and b-B n P n structures are in the range from À4.723 to À4.803 eV per atom, which is much better or comparable to that of other 2D materials such as phosphorene, 68 carbon phosphide, 69 Al 2 C sheet, 70 and silicene, 71 having 3.48, 5.32, 3.94, and 3.55 eV per atom, respectively, indicating the very good stability of these materials. Furthermore, the binding energy of the predicted clusters was compared with that of different 2D-BP allotropes. The calculated binding energies for the a-B n P n and b-B n P n structures is comparable with the binding energy of 2D a 0 -BP (4.91 eV per atom), a 1 -BP (4.82 eV per atom), b 0 -BP (4.72 eV per atom), b 1 -BP (4.54 eV per atom), and g-BP (4.44 eV per atom). 72 Moreover, the binding energies for these new congurations are even close to the already synthesized B 40 , 73 M@B n , 66 and B 24 clusters. 64 Hence, it can be concluded that the a-B n P n and b-B n P n structures have excellent stability and possible experimental synthesis. As seen in Table 2, the a-B n P n structures are more stable than the corresponding b-B n P n structure and the maximum binding energy of À4.803 eV per atom belongs to the a-B 24 P 24 structure.

Dynamic and chemical stability
Next, the dynamic stability of the predicted clusters was also investigated. The vibrational frequencies at the gamma points for the most stable structures as well as the rst low-lying isomers were calculated. Based on these calculations, no negative vibration frequency could be found, conrming that the most stable and the rst low-lying structures are dynamically stable. Furthermore, the thermal stability of the most stable structure and the rst low-lying energy isomers of the studied B n P n clusters was performed using ab initio molecular dynamic simulations (AIMD) at three different (300, 600, and 900 K) temperatures by the NVT-ensemble with the Nosé-Hoover thermostat. As represented in Fig. 3 and 4, the predicted structures equilibrated at 300 K from the AIMD simulations. Furthermore, there is no broken bond and the initial clusters are well maintained during the AIMD simulation for 5 ps with a time-step of 1 fs. As shown in these gures, the total potential energy only uctuates around a certain constant magnitude, which indicates the most stable structures of the B n P n clusters and the rst low-lying isomers are thermodynamically stable and could hold their structural integrity at room temperature. It can be seen from Fig. 3 and 4 that aer heating both a-B n P n and b-B n P n structures, the basic structure of the clusters is also well maintained at 600 and 900 K. Furthermore, the thermal stability of the a-B 24 P 24 cluster at a higher temperature such as 1200 K was also examined (Fig. S4 †). Based on the AIMD simulations, when the temperature reaches as high as 900 K, the bond elongation of the 6-membered ring occurs. It can be anticipated that as the temperature increases, bond stretching is also increased. Hence, the upper limit of temperature that the structures can tolerate is 900 K and the bond stretching occurs between 900 and 1200 K.
Chemical stability in air and other environments is an important issue that may limit the future applications of materials. Furthermore, a major problem with atoms located on the surface of the materials is their reaction with oxygen, nitrogen, and water in air; hence, the surface chemistry of these materials is essential in determining the intrinsic properties. 74,75 Some phosphorus-based materials such as 2Dphosphorene and black phosphorus are known to oxidize and break down in ambient conditions because of the lone atomic pairs of the P atoms at the surface. [76][77][78][79][80] Generally, P atoms favor the formation of the trigonal pyramidal sp 3 bonded conguration. 81 It is found that in the a-B 12 P 12 cluster, each P atom with sp 3 hybridization bonding with two P and one B neighbor atoms has one lone pair of electrons; hence, it motivates us to investigate the effects of oxygen, nitrogen, and water-gas molecules on the stability of the a-B n P n structures.
The AIMD simulation was used to investigate the interaction between a-B n P n and gaseous phase O 2 , N 2 , and H 2 O molecules at room temperature. As an initial model, randomly, eight O 2 , N 2 , and H 2 O molecules were added around the a-B 12 P 12 anda-B 24 P 24 clusters with 1 Â 1 Â 1 unit cells at an initial distance of 3.5-4 A from the cluster surface (see Fig. 5 and 6). The lattice parameters of the optimized a-B n P n structures at the PBE + D3 level are a ¼ 19.4 A, b ¼ 19.4 A, and c ¼ 19.4 A. The AIMD simulations were performed for 7 ps with a time step of 1.0 fs under 300 K. Based on our AIMD simulations, it can be seen that generally the a-B n P n clusters have relatively high stability in the presence of O 2 , N 2 , and H 2 O molecules. Obviously, aer 7 ps of contact, some O 2 molecules dissociated into O atoms and chemisorbed onto the a-B 12 P 12 surface. The P-O bonds with a bond length of about 1.524 A are short and could be polar due to the difference in the electronegativity between the P and O atoms. Apart from the negative outcomes due to the a-B 12 P 12 interaction with the O 2 molecule, controlled passivation might lead to new stable structures for gas sensing. It should be reminded that the basic structure of the a-B 12 P 12 cluster is also well maintained aer 7 ps of contact. In contrast, the phenomena of O 2 dissociation has not been found on the a-B 24 P 24 surface in the AIMD simulation (see Fig. 6), which indicated the high oxidation resistance of the a-B 24 P 24 structure. This result is very exciting and is different from the result we expected. Furthermore, the simulations identify that no interaction is observed between the a-B 12 P 12 /a-B 24 P 24 structures with N 2 and H 2 O molecules, which indicates the high stability of this type of structure. To conrm the chemical stability of the a-B n P n clusters, the adsorption of O 2 , N 2 , and H 2 O gas molecules onto the a-B n P n surfaces was performed. The adsorption energy (E ads ; eV per gas molecule) can be calculated as follows.
E ads ¼ (E cluster-gas À (E cluster + nE gas )/n where n is the number of gas molecules per unit cell. The E cluster-gas term stands for the total energy of the a-B n P n clusters aer gas adsorption. The E cluster and E gas terms are the total energies of the isolated a-B n P n cluster and the gas molecule, respectively. According to the above expression, the negative E ads exhibits that complex formation is thermodynamically Table 2 Symmetries, calculated binding energy per atom (E bin , eV per atom), HOMO-LUMO energy gaps (D g , eV), and VBM and CBM (eV) of the (BP) n clusters. The structural and energetic properties are calculated using the PBE-D3 method, while the HSE06 functional was used to calculate the electronic properties (D g ), VBM, CBM, and summation Bader partial charge on all B and P atoms (esu) of the HSE06 functional based on the PBE-D3 relaxed structures favorable. The calculated adsorption energy of the studied complexes is plotted in Fig. S5. † As seen from this gure, the adsorption energy of the O 2 molecules onto the a-B 12 P 12 surface is favorable (À42.93 kJ mol À1 per O 2 molecule). According to the adsorption energy values, it is found that O 2 physisorbs over a-B 12 P 12 with an adsorption energy of about À43 kJ mol À1 . These results are in good accordance with the results of the AIMD simulation, which were used to investigate the chemical stability a-B n P n and O 2 , N 2 , and H 2 O gas molecules. Hence, it can be concluded that the structural integrity of a-B n P n can be retained perfectly in the environment of gaseous phase O 2 , N 2 , and H 2 O at room temperature.

Electronic and optical properties
Based on the relaxed geometries, the electronic properties of the most stable structure, along with the rst low-lying isomer of the B n P n structures, were also investigated. The HOMO-LUMO (H-L) energy gap (E gap ) is a fruitful parameter for examining the kinetic stability in molecules. The E gap is the energy needed to excite the electron from the HOMO to the LUMO, where a large Fig. 3 The snapshots of the most stable structures for a time of 5 ps of (a) a-B 12 P 12 and (b) b-B 12 P 12 . The temperature of the AIMD was set to 300, 600, and 900 K.
energy gap implies more kinetic stability and less reactivity. The HSE06 energy gaps of the ground-state structures and the rst low-lying isomer of B n P n are summarized in Table 2. As represented in Table 2, the H-L energy gaps of the a-B n P n and b-B n P n clusters are in the range of $0.67-3.2 eV, while the energy gaps in the b-B n P n fullerene-like structures are higher than those of the corresponding a-B n P n structures. It is revealed that a-B 12 P 12 and a-B 24 P 24 are semiconductor structures with the HSE06 band gap of 0.67 and 0.62 eV respectively, which is comparable with the band gap of the 2D-BP structure. 28 In contrast, b-B 12 P 12 and b-B 24 P 24 fullerene-like structures have large band gaps of 3.2 and 2.98 eV, respectively. The small band gap value suggests Paper Nanoscale Advances modest chemical stability while the sizeable HOMO-LUMO energy gap is representative of the sufficient chemical stability of the studied clusters. The charge density distribution of the HOMO and LUMO for both the a-B n P n and b-B n P n structures are presented in Fig. 7. As seen, the HOMO and LUMO orbitals are generally composed of p-orbitals of phosphorus and boron atoms, and they have been distributed over the boron units and B-P bonds. In case of a-B n P n , each P atom with sp 3 conguration has one lone pair of electrons in the fourth orbital (bonded with three neighbors of two P and one B atoms). 81 Phosphorus has a higher electronegativity than the boron atom; hence, it is expected that the charge will be transferred from the B to the P atoms in the BP clusters. Bader charge analysis method is evaluated using the HSE06 functional 82 and the corresponding values are collected in Table 2. As expected, boron atoms become slightly positively charged with the loss of electron. The amount of total charge transfer from the B atoms to the P atoms in the b-B n P n structures is more signicant than the corresponding a-B n P n structures. Due to the charge transfer between B and P atoms, it is demonstrated that the B n P n clusters have polar properties, especially for the b-B n P n structures. It can be concluded that charge transfer between the B and P atoms demonstrates that electrostatic interactions have a signicant effect on the stability of these types of structures. To understand more about the nature of bonds in the B n P n structures, the electron localization function (ELF) calculated by VASP is also investigated. ELF can represent electron localization in a molecule and in the solid-state, and its values span from 0 to 1, which illustrates the degree of electron localization. Generally, the considerable value ($1) of ELF is consistent with the strong covalent bonds or lone pair of electrons. On the other hand, a small value of ELF points out ionic or metallic bonds, and it shows the low electron density localization. Fig. 8 and S5 † represent the ELF (isovalue ¼ 0.8 au) of the B n P n structures. From these results, it can be realized that electron localization around the P atoms is slightly more than that of the B atoms, which in good agreement with Bader charge analysis. Moreover, the line plots of the ELF values for the B-P bonds indicates the covalent nature of these bonds in the B n P n structures (Fig. S6 †).
The high stability and moderate band gap feature of the investigated B n P n clusters motivate us to explore the optical properties of the entitled structures. The optical properties of the B n P n clusters have also been studied using the HSE06 functional. To investigate the performance under light, the optical absorption coefficient a was calculated according to the equation below. 75 where 3 1 (u) and 3 2 (u) refer to the real and imaginary parts of the complex dielectric functions, respectively. As displayed in Fig. 9, B n P n materials exhibit notable absorption in both visible and ultraviolet range with the absorption coefficient larger than 10 5 cm À1 , compared to many 2D-materials used as electronic Paper Nanoscale Advances and optoelectronic materials. 84 According to Fig. 9, the optical absorption spectrum of both the structures illustrates a higher degree of absorption in the UV and visible region, which suggests the ability of the a-B n P n and b-B n P n structures to capture ultraviolet and visible light harvesting. In comparison, a-B n P n exhibits stronger optical absorption than those corresponding to the b-B n P n structures in the visible region. Besides, we focused on the potential applications of these compounds for photocatalytic water-splitting. First, the HSE06 position of the valence band maximum (VBM) and the conduction band minimum (CBM) was calculated from the difference between the electrostatic potential at the vacuum region and is shown in Fig. 10. As discussed above, the b-B 12 P 12 and b-B 24 P 24 structures have a band gap of 3.2 and 2.98 eV, respectively, which is located in the visible-light region. The optical absorption of the b-B 12 P 12 and b-B 24 P 24 fullerene-like structures, as shown in Fig. 9, represents that the b-B 12 P 12 and b-B 24 P 24 clusters have expectantly optical absorption in the   Optical absorption spectrum of (BP) n with n ¼ 12 and 24 clusters calculated by the HSE06 functional. The energy range corresponding to visible light is indicated by vertical red-dashed lines. B 12 P 12 and b-B 24 P 24 clusters meet both conditions for photocatalytic water splitting. As seen from Fig. 10, their valence bands lie at more positive potentials than the water oxidation potential (O 2 /H 2 O potential) and their conduction bands (CBMs) are more negative than the hydrogen reduction potential (H + /H 2 potential). These results clearly show that the b-B n P n structures have great potential for applications as a visible-light photocatalyst for water splitting. In contrast, the strong absorption coefficient of the a-B 24 P 24 (which can reach a higher order of 5 Â 10 5 cm À1 ) and a-B 12 P 12 systems is comparable to that of the organic perovskite solar cell. These results reveal that a-B 24 P 24 and a-B 12 P 12 have great potential for applications as optical absorbent materials in solar cells and optoelectronic devices.

NLO properties
Based on the research conducted to investigate the ability of (XY) n fullerene-like clusters as NLO materials, 31,38,85,86 it has motivated us to investigate the second-order NLO properties of both a-B n P n and b-B n P n structures. Second-order NLO materials that can alter the wavelength of lasers are of considerable interest because of the potential of utilization of such materials in optical switching, optical computing, and photonic technology. 31,38,85,86 Quartz crystal (a-SiO 2 ) is a rst second-order NLO material, as demonstrated by Franken et al. in 1961. 87,88 In most recent studies, the ability of metal encapsulated (XY) n nano-cages (X ¼ B, Al, and Y ¼ N, P) as NLO materials has been considered. They found that the (XY) n nano-cages have zero dipole moment and zero hyperpolarizability. It should be noted that none of these studies have examined the second-order NLO properties. Herein, for both a-B n P n and b-B n P n structures, the average polarizability (a), rst hyper-polarizabilities (b 0 ), and second hyperpolarizability (g) were calculated at the CAM-B3LYP/def2-tzvp level of theory using the Gaussian 09 package. 89 As seen in Table 3, the mean polarizabilities (a) of a-B n P n are higher than those of the corresponding b-B n P n structures. The maximum value of a is related to the a-B 24 P 24 structure with 965.24 a.u., which is higher than that of the many encapsulated metal-(XY) n nano-cages. 31,86 Furthermore, our calculations proved that both a-B n P n and b-B n P n structures have nearly zero rst hyper-polarizabilities (b 0 ). The second hyperpolarizabilities (g) of a-B n P n and b-B n P n are also summarized in Table 3. The second hyper-polarizabilities of all the predicted structures are higher than the experimental values of g for C 60 and C 70 fullerenes, 90,91 and comparable to that of the best NLO materials reported in the literature. On the other hand, a-B 12 P 12 and a-B 24 P 24 have larger g values than the corresponding b-B n P n clusters. The g value of a-B 24 P 24 is 29 times greater than the second hyper-polarizability value of the b-B 24 P 24 clusters. In conclusion, the strong interactions, effective intermolecular charge transfer, and cluster shape affect the magnitude of the second hyper-polarizabilities of the a-B 12 P 12 and a-B 24 P 24 structures. Furthermore, the a-B 12 P 12 and a-B 24 P 24 clusters show not only high stability but also have large second hyperpolarizability. Hence, they are anticipated to be excellent potential unprecedented candidates for second-order NLO materials.

Conclusion
In this work, we explored the new stable structure of B n P n (n ¼ 12, 24) clusters with a modern method of crystal structure prediction, USPEX conjugated with spin-polarized rstprinciples calculations. Our prediction demonstrates that the most stable structures of B 12 P 12 and B 24 P 24 are a-B 12 P 12 and a-  B 24 P 24 . Based on the relaxed geometries of the a-B n P n structures, each phosphorus atom is doped into a boron atom while the B atoms form a B n unit. In contrast, the rst low-lying b-B n P n isomers of the (BP) n clusters are shaped like fullerenes, where each boron atom is bonded to each phosphorus atom to make cage structures. The calculated binding energies for these novel structures are in the range from À4.723 to À4.803 eV per atom, which is much better or comparable to that of other 2D materials and clusters, indicating their excellent stability and possible experimental synthesis. Furthermore, a-B n P n structures are more stable than the corresponding b-B n P n structures, and the maximum binding energy of À4.803 eV per atom belongs to the a-B 24 P 24 structure. Based on our AIMD simulations, it was found that a-B n P n structures have relatively high stability in the presence of O 2 , N 2 , and H 2 O molecules, in particular, a-B 24 P 24 shows superior oxidation resistance. In the case of thermal stability, a-B n P n and b-B n P n exhibit better thermal stability, where the upper limit temperature that the structures can tolerate is 900 K. It is revealed that a-B 12 P 12 and a-B 24 P 24 are semiconductor structures with the HSE06 band gap of 0.67 and 0.62 eV, respectively, which is comparable with the band gap of the 2D-BP structure. In contrast, b-B 12 P 12 and b-B 24 P 24 fullerene-like structures have large band gaps of 3.2 and 2.98 eV, respectively. The calculated Bader charges illustrate their ionic characters with charge transfers from the B to P atoms. The electronic properties of these novel compounds illustrate a higher degree of absorption in the UV and visible-region with an absorption coefficient larger than 10 5 cm À1 , which suggests a wide range of opportunities for advanced optoelectronic applications. The b-B n P n phase has suitable band alignments in the visible-light excitation region, which will produce enhanced photocatalytic water splitting. On the other hand, a-B n P n structures with high absorption coefficient exhibit large second hyperpolarizability, which are expected to have excellent potential as second-order NLO materials.

Conflicts of interest
There are no conicts to declare.