# Discovery of Novel Silicon Allotropes with Optimized Band Gaps to Enhance
Solar Cell Efficiency through Evolutionary Algorithms and Machine Learning

Mostafa Yaghoobi  Mojtaba Alaei  Mahtab Shirazi  Nafise Rezaei  Stefano de
Gironcoli  Department of Physics, Isfahan University of Technology, Isfahan
84156-83111, Iran  Scuola Internazionale Superiore di Studi Avanzati, Trieste,
Italy  CNR-IOM DEMOCRITOS, Istituto Officina dei Materiali, Trieste, Italy
Skolkovo Institute of Science and Technology, 121205, Bolshoy Boulevard 30,
bld. 1, Moscow, Russia

###### Abstract

In the pursuit of advancing solar energy technologies, this study presents 20
direct and quasi-direct band gap silicon crystalline semiconductors that
satisfy the Shockley-Queisser limit, a benchmark for solar cell efficiency.
Employing two evolutionary algorithm-based searches, we optimize structures
and calculate fitness function using the DFTB method and Gaussian
approximation potential. Following the preselection of structures based on
energy considerations, we further optimize them using PBEsol DFT.
Subsequently, we screen the structures based on their band gap, employing a
DFTB method tailored for band gap calculation of silicon crystals. To validate
the structural stability, we employ phonon analysis via linear regression
algorithm applied to PBEsol DFT data. To ensure accurate band gap
determination, we employ HSE and GW methods. Significantly, the structures
unveiled in this study are of great importance due to their proven stability
from both mechanical and dynamic perspectives. Furthermore, the ductility and
low density of certain structures enhance their potential application.
Moreover, the study harnesses various machine learning algorithms to develop a
predictive model for the band gap energy of these silicon structures. Input
data for machine learning models are derived from structural MBTR and SOAP
descriptors, as well as DFT outputs. Notably, the results reveal that features
extracted from DFT outperform the MBTR and SOAP descriptors.

##  1 Introduction

In 2016, the Paris Agreement established a goal of limiting the global average
temperature increase to well below $2^{\circ}$C above pre-industrial levels,
with a further aim of restricting it to $1.5^{\circ}$C. Fulfilling this
crucial mission requires the development of alternative, clean energy sources
that can reduce our dependence on fossil fuels. Solar energy is a prime
candidate for a sustainable and eco-friendly energy source, and photovoltaic
(PV) cells are key to converting solar energy into electricity. However, the
efficiency of PV cells is limited by the properties of the materials used to
construct them, particularly their band gap. Moreover, employing direct and
quasi-direct bandgap semiconductors in solar cells results in a diminished
thickness of the solar cell. In these materials, there is enhanced absorption
of incident photons, leading to increased cell efficiency [1]. According to
the Shockley–Queisser limit [2], solar cell efficiency surpasses 30 percent
when the band gap of the material used in the cell falls within the range of
0.93 to 1.61eV, with the highest efficiency occurring at the band gap of 1.34
eV. While silicon has historically been the primary material for PV cells due
to its abundance and stability, its indirect band gap energy of 1.1 eV poses
limitations on its efficiency. Therefore, finding silicon structures that have
a direct and quasi-direct band gap energy and are within the Shockley–Queisser
limit will be an important step in the development of solar cells and in
increasing their efficiency.

Predicting crystal structures has always been a significant challenge in
material science. In Ref. 3, Fan et al. proposed two phases of metastable
silicon structures by substituting the Si atom in the carbon structures. Zhang
et al. [4] proposed four semiconducting silicon allotropes with direct and
quasi-direct band gaps between 1.193 and 1.473 eV, characterized by low
energy. Chun-Xiang Zhao et al. [5] have introduced a new silicon allotrope
with a direct band gap of 0.591eV. Meanwhile, Wang et al. [6] suggested six
metastable silicon structures using the particle swarm optimization algorithm.
Additionally, using the inverse band structure design approach based on the
particle swarm optimization algorithm, H. J. Xiang et al. were able to predict
a 20-atom silicon phase with a quasi-direct gap of 1.55 eV [7]. Moreover, in a
study by Q. Wei et al. [8], researchers explored crystalline silicon
structures with direct gaps by replacing carbon atoms with silicon in the
existing carbon phases, yielding six silicon allotropes with direct band gaps.
Also, in Refs. 9 and 10, researchers introduced several metastable silicon
crystal structures with direct and quasi-direct band gaps within the range of
1 to 1.8 eV, employing the ab initio minima hopping method.

In this study, we used the power of the evolutionary algorithm (EA) [11, 12,
13] to predict novel phases of silicon crystal structures, a method that has
recently captured the attention of researchers. Through two distinct EA-based
searches, we explored silicon structures ranging in size from 2 to 45 atoms.
In the first search, we employed the self-consistent-charge density-functional
tight-binding (DFTB) [14, 15] method to calculate the fitness function, which
in this study is the enthalpy of the structures. This choice aimed to ensure
the stability of the predicted structures. In the second search, the Gaussian
approximation potential (GAP) [16] was used to calculate the fitness function.
The GAP potential significantly accelerated the enthalpy predictions, thereby
facilitating the exploration of larger structures.

Subsequently, our examination focused on identifying silicon structures that
align with the Shockley–Queisser limit, thereby demonstrating the potential
for enhanced solar cell efficiency. Through a comprehensive analysis of the
dynamical and mechanical stability of these structures, we sought to guarantee
their practicality and reliability. Ultimately, we meticulously screened and
identified a collection of 20 unique and promising silicon structures, each
holding the potential to reshape the landscape of solar cell technology.

The paper is structured as follows: we begin by detailing our methodology,
casting the foundation for the subsequent results. We then unveil the silicon
structures we have discovered and describe their mechanical and electronic
properties. Subsequently, our discussion addresses the implications of these
findings and highlights the most promising structures for solar cell
applications. Finally, we delve into the prediction of band gap energy for
these structures, utilizing various machine learning(ML) models.

##  2 Computational Details

In the pursuit of crystal structure prediction, we harnessed the well-
established USPEX code [17], which has demonstrated promising results and is
extensively employed in materials science researches [18, 19].

As mentioned, Our study encompassed two distinct searches, each wielding
distinct computational approaches. In the first search, to calculate the
enthalpy of the structures we used DFTB with pbc Slater-Koster [20] (SK),
executed through the DFTB+ software package [21], which offers a
computationally efficient alternative to first-principles DFT methods. In the
subsequent search, we replaced DFTB with GAP [22]. GAP’s formalism rests on
the Gaussian process regression (GPR) and smooth overlap atomic position
(SOAP) [23] kernel. To implement the GAP and incorporate it with USPEX, we
used the LAMMPS code [16, 24].

In the first search, we set the initial population size to $250$, the
population of each generation to $200$ structures, the number of generations
to $100$, and a stopping criterion, whereby the algorithm would terminate if
there was no improvement in the best structures after 30 generations.
Furthermore, in our pursuit of more optimized structures, we employed the
USPEX once more to predict atomic structures comprising 2 to 20 atoms. In this
time, we employed the structures procured during the preceding search as
initial seeds within the genetic algorithm framework. We carried out two
generations of the genetic algorithm, with a population size of $50$ in each
generation. For fitness function calculations and structural relaxation, we
employed DFT via Quantum Espresso (QE) [25].

In the second search, we set the initial population size to 300 and the
population of each generation to 250 structures. We specified the number of
generations to 200, and a stopping criterion, whereby the algorithm would
terminate if there was no improvement in the best structures after 30
generations. Overall, according to the enthalpy and volume of the structures,
the first search led to the production of 536 structures, and the second
search led to the production of 827 structures.

To ensure that equilibrium positions were reached for each structure, all
structures (1263 Si structure) were relaxed by QE using projector augmented
wave method (PAW) with higher accuracy. We set PBEsol functional for exchange-
correlation appximation, a uniform reciprocal space mesh point spacing of 0.15
(1/Å), 40 Ry kinetic energy cutoff for wavefunctions, and 400 Ry kinetic
energy cutoff for charge density.

After relaxing the structures, we investigated the band gap energy of the
structures to find suitable structures for solar cells. Despite the
satisfactory accuracy of DFT in material energy calculations, this method
leads to underestimation in band gap energy prediction. In order to overcome
the band gap problem, methods such as GW [26] and the Heyd-Scuseria-Ernzerhof
(HSE) hybrid functional [27] have been suggested, although they are
computationally demanding. More recently, researchers have turned to the DFTB
method, leveraging the Hamiltonian and overlap matrix parameters of atomic
orbitals, which are pre-tabulated and calculated using empirical methods. This
approach offers a faster and more cost-effective means to obtain the band
structure of materials. In our study, we utilized the DFTB method, employing
si-band Slater-Koster [28] fitted for band structure calculations, to compute
the band structure of the crystal structures. To assess the accuracy of the
DFTB method in comparison to the GW and HSE hybrid functional methods, we
selected four structures previously identified as promising candidates for use
in solar cells and computed their band structure within the self-consistent GW
using the Yambo code [29]. The resulting band gap energy values for each
structure and each method (HSE, GW and DFTB) are presented in the Table 1.

Table 1: Values of the band gaps in eV that obtaind by HSE, GW and DFTB. The
structures are extracted from Ref [41].

structure |  HSE[41] | GW | DFTB
---|---|---|---
mC12 | 1.24 | 1.33 | 1.31
oC12 | 0.96 | 1.03 | 1.19
tl16 | 1.25 | 1.33 | 1.09
tp16 | 0.91 | 0.90 | 1.19

The table clearly illustrates that the DFTB method based on si-band SK
exhibits a good level of accuracy when compared to the other two methods.
Importantly, DFTB stands out for its significantly reduced computational time
requirements. It is worth mentioning that we use Wannier90 [30] to obtain the
HSE band structures from HSE scf calculations of QE.

In our research, we also conducted a comprehensive assessment of the
mechanical and dynamic stability of the screened crystal structures. To
achieve this, we employed DFT calculations using the THERMO-PW package code
[31]. The DFT calculations were performed, using the PBEsol pseudopotential
[32]. Dynamically stable structures can withstand small atomic vibrations, a
property we examined by studying their phonon properties. We utilized the ML-
based hiphive code [33] to check the phononic properties of the structures.

Furthermore, We conducted a study aimed at developing an ML model to predict
the band gap energy of silicon crystal structures. In our investigation, we
explored various ML algorithms to determine their effectiveness in predicting
the band gap of these structures. The algorithms considered in our analysis
included Multilayer Perceptron Neural Network that we simply denote by NN,
Kernel Ridge Regression (KRR), Support Vector Regression (SVR), Convolutional
Neural Networks (CNN), as well as Decision Tree (DT), Random Forest (RF), and
XGbosst. We implemented NN and CNN using the TensorFlow framework [34], while
for SVR, DT, RF, and XGBoost, we utilized the Scikit-Learn library [35] in
Python. In this study, we use two important structural descriptors, Many-body
tensor representation(MBTR) [36] and SOAP [37], to construct suitable features
for ML models. To generate these descriptors for the crystal structures, we
utilized the Dscribe computational package [38]

##  3 Results

###  3.1 Structure Screening

In our initial search, we identified 148 semiconductor structures featuring
band gaps ranging from 0.0949 to 2.014 eV. In the second search, it turns out
that there are 256 semiconductor structures with band gap between 0.007 and
2.183eV. To ensure the uniqueness and distinctiveness of our selection, we
conducted a thorough examination of the structural geometry. This scrutiny led
to the elimination of similar structures and those that proved identical to
the silicon diamond structure. This process was done by comparing the Valle-
Oganov fingerprint [39] of each structure with other structures and the
silicon diamond structure. Therefore, the number of unique semiconductor
structures was reduced to $88$ in the first search and to 108 structures in
the second search. The band gap histogram of the structures for both searches
depicted in Figure 1. The average band gap energy value in the initial search
equals 1.06 eV, while in the subsequent search, it is 0.72 eV.

Figure 1: The histogram of the band gap energies for the structures in both
searches.

Further narrowing our focus, we identified 85 structures that have band gaps
between 0.93 and 1.61 eV. Screening these structures for those whose band gap
is direct and quasi-direct, we found a subset of 20 structures suitable for
further consideration.

Figure 2: The ball-and-stick representations of the candidate crystal
structures. The space group of these structures is presented in Table 4.

In addition, we conducted a geometric analysis of all 20 screened final
structures through Spglib[40]. By increasing the tolerance for the Cartesian
coordinate distances between atoms, we achieved geometric refinement,
resulting in symmetric structures. Following this refinement, we performed
another optimization of the structures using PBEsol DFT. The ball-and-stick
drawings of the crystal structures are presented in Figure 2.

Furthermore, an examination of the fingerprint of all 20 structures was
conducted in comparison to the pure silicon structures introduced in the
Material Project database. This analysis revealed that Si16(II) had previously
been documented on the Material Project platform. Moreover, structures
Si12(II) and Si16(III) exhibited a noteworthy similarity to structures oC12
and tl16, respectively, as previously elucidated in Ref. 41.

It is worth highlighting that among these structures, Si40, Si12(II),
Si16(III), Si20(III), Si16(II), Si16(I), Si36(II), and Si36(I) hold particular
significance due to their classification as porous crystals. These structures
exhibit void tunnels, a distinct feature in their composition. The void
tunnels in Si40, Si12(II), Si36(II), Si16(I), Si16(II), and Si20(III)
structures are formed by connecting 4, 5, and 6-membered rings, and are
observable from one side. However, in the case of Si36(I), the void tunnels
exist on two sides, and they are formed by connecting 4-membered rings.

Direct and quasi-direct semiconductors are particularly advantageous in the
fabrication of thin-film solar cells, obviating the need for a thick layer
required in the use of indirect band gap semiconductors [1]. To evaluate the
band gap properties of various semiconductor structures, we conducted an in-
depth analysis of their band structures. The HSE and GW band structures can be
found in Supplementary Information (SI). The results of this analysis are
summarized in Table 2, wherein DFTB, HSE, and GW results were showed for a
comprehensive examination.

Table 2: The band gap energies of the structures determined via DFTB, HSE, and GW calculations. Additionally, the band gap energies of the structures predicted using ML models including CNN, XGBoost, and RF, in combination with DFT outputs features. In the second column, we indicate the structure type by specifying whether it is a direct(D) or quasi-direct(QD) semiconductor. structures | type | band gap | ML band gap
---|---|---|---
DFTB | HSE | GW | CNN | XGBoost | RF
Si12(I) | QD | 1.28 | 1.46 | 1.08 | 1.08 | 1.27 | 1.17
Si12(II) | D | 1.24 | 1.16 | 0.75 | 1.24 | 1.21 | 1.15
Si12(III) | QD | 1.49 | 1.63 | 1.34 | 1.56 | 1.49 | 1.47
Si14 | D | 1.13 | 1.06 | 0.85 | 0.91 | 1.11 | 1.01
Si16(I) | D | 0.97 | 1.15 | 0.93 | 1.11 | 0.90 | 0.90
Si16(II) | QD | 0.95 | 0.92 | 0.73 | 0.72 | 0.68 | 0.59
Si16(III) | D | 1.16 | 1.35 | 1.19 | 1.24 | 1.11 | 1.17
Si16(IV) | D | 1.18 | 1.31 | 1.11 | 1.11 | 1.07 | 1.08
Si18(I) | QD | 1.43 | 1.6 | 1.42 | 1.37 | 1.27 | 1.34
Si18(II) | D | 1.03 | 1.46 | 1.30 | 1.24 | 1.05 | 1.20
Si20(I) | D | 1.22 | 1.31 | 1.03 | 1.27 | 1.21 | 1.21
Si20(II) | D | 1.26 | 1.45 | 1.20 | 1.31 | 1.25 | 1.23
Si20(III) | QD | 1.16 | 1.31 | 1.06 | 1.20 | 1.21 | 1.16
Si22 | D | 1.4 | 1.57 | 1.34 | 1.49 | 1.35 | 1.24
Si28 | QD | 1.38 | 1.50 | 1.23 | 1.31 | 1.36 | 1.29
Si36(I) | QD | 1.15 | 1.43 | 1.28 | 1.47 | 1.25 | 1.27
Si36(II) | D | 1.14 | 1.28 | 1.17 | 1.10 | 1.05 | 1.08
Si36(III) | D | 1.13 | 1.34 | 1.20 | 1.17 | 1.12 | 1.14
Si40 | QD | 1.34 | 1.35 | 1.35 | 1.22 | 1.20 | 1.24
Si44 | QD | 1.38 | 1.24 | 1.07 | 1.32 | 1.36 | 1.28
Si${}_{\text{diamond}}$ | ID | 1.14 | 1.28 | 1.10 | 1.11 | 1.26 | 1.23

The DFTB has a 14% error, and the HSE has a 20% error when compared to GW
results. This suggests that the DFTB method is reliable for estimating the
band gap. In addition, Table 2 contains the values of band gap energy
predicted by three of the best ML models, including (CNN, XGBoost, and RF)+DFT
output. According to our findings, 11 structures including Si16(III),
Si12(II), Si16(I), Si22, Si36(II), Si14, Si18(II), Si36(III), Si20(II),
Si20(I), and Si16(IV), exhibit direct band gaps, while the remaining
structures fall into the category of quasi-direct semiconductors. To ascertain
whether a structure qualifies as a quasi-direct semiconductor, we compute the
direct band gap (E${}_{g}^{d}$) and the band gap (Eg) If the result of
E${}_{g}^{d}$ \- Eg is less than 0.15 eV, the system is classified as a quasi-
direct band gap semiconductor [42]. The corresponding values of Eg,
E${}_{g}^{d}$, and E${}_{g}^{d}$ \- Eg for the quasi-direct semiconductors
obtained through HSE method are presented in the Table 3.

Table 3: Band gap energy (Eg) and direct band gap (E${}_{g}^{d}$) values of the structures that were known quasi-direct semiconductor based on HSE calculations. structure | E${}^{d}_{g}$ | Eg | E${}^{d}_{g}$ \- Eg
---|---|---|---
Si12(I) | 1.50 | 1.46 | 0.043
Si12(III) | 1.64 | 1.63 | 0.008
Si16(II) | 1.06 | 0.92 | 0.139
Si18(I) | 1.65 | 1.60 | 0.044
Si20(III) | 1.32 | 1.31 | 0.005
Si28 | 1.50 | 1.50 | 0.003
Si36(I) | 1.52 | 1.43 | 0.094
Si40 | 1.39 | 1.35 | 0.043
Si44 | 1.24 | 1.23 | 0.006

It is important to emphasize that there are cases where the results of GW and
HSE analysis differ. Specifically, an examination of the GW calculations for
the Si16(IV) structure reveals it to be a quasi-direct band gap with a direct
gap of 1.13 eV. However, HSE results designate this structure as a direct
semiconductor. In the case of Si18(II), HSE calculations categorize it as a
direct semiconductor, while GW analysis identifies it as a quasi-direct
semiconductor with a direct gap of 1.37 eV. A similar pattern is observed in
the assessment of Si22, where the direct gap is determined to be 1.36 eV
through GW calculations. Furthermore, Si28 and Si20(III) structures are
identified as direct semiconductors according to GW, but HSE calculations
classify them as quasi-direct semiconductors.

###  3.2 Structural Properties

Table 4 presents the data encompassing bulk modulus ($B$), shear modulus
($G$), Young modulus ($E$), Poisson ratio ($n$), and the ratio of bulk modulus
to shear modulus($B$/$G$) in the Voigt-Reuss-Hill average approximation.

Table 4: The table comprises bulk modulus ($B$), Young modulus ($E$), shear
modulus ($G$) measured in gigapascals (GPa), Poisson ratio ($n$), the ratio of
bulk to shear modulus ($B$/$G$), density ($\rho$) expressed in grams per cubic
centimeter (g/cm3), and relative enthalpy ($\Delta H$) per atom expressed in
electronvolts (eV). Last column shows the space group of structures.

structure |  $B$(GPa) | E(GPa) | G(GPa) | $n$ |  $B$/$G$ | $\rho(g/cm^{3})$ | $\Delta H\text{(eV per atoms)}$ | space group
---|---|---|---|---|---|---|---|---
Si12(I) | 68 | 82 | 32 | 0.297 | 2.142 | 2.1428 | 0.2165 | C2/m
Si12(II) | 83 | 132 | 54 | 0.233 | 1.54 | 2.1650 | 0.0924 |  C2221
Si12(III) | 87 | 134 | 54 | 0.243 | 1.617 | 2.2866 | 0.0028 |  Pmn21
Si14 | 95 | 154 | 63 | 0.229 | 1.514 | 2.5256 | 0.0104 | C222
Si16(I) | 89 | 131 | 52 | 0.256 | 1.720 | 2.2727 | 0.0644 | Pmmm
Si16(II) | 85 | 132 | 53 | 0.241 | 1.595 | 2.2452 | 0.0710 | C2/m
Si16(III) | 80 | 108 | 42 | 0.272 | 1.881 | 2.0726 | 0.1419 | I4/mcm
Si16(IV) | 85 | 115 | 45 | 0.276 | 1.900 | 2.2054 | 0.0083 | Cmcm
Si18(I) | 87 | 137 | 55 | 0.237 | 1.571 | 2.3230 | 0.0691 | P-1
Si18(II) | 76 | 112 | 45 | 0.254 | 1.704 | 2.2050 | 0.0145 | P-1
Si20(I) | 81 | 122 | 49 | 0.247 | 1.648 | 2.3782 | 0.0121 |  P43212
Si20(II) | 87 | 134 | 54 | 0.244 | 1.621 | 2.2285 | 0.0039 | Imm2
Si20(III) | 85 | 134 | 54 | 0.236 | 1.562 | 2.2242 | 0.0707 | C2/m
Si22 | 83 | 115 | 45 | 0.268 | 1.829 | 2.1605 | 0.1382 |  P42/mcm
Si28 | 87 | 132 | 53 | 0.246 | 1.641 | 2.2896 | 0.0045 | C2/m
Si36(I) | 51 | 51 | 18 | 0.334 | 2.675 | 1.6026 | 0.3854 | R3m
Si36(II) | 83 | 132 | 53 | 0.235 | 1.560 | 2.2913 | 0.0876 | C2/m
Si36(III) | 82 | 117 | 46 | 0.263 | 1.777 | 2.1329 | 0.0091 | Ibam
Si40 | 51 | 62 | 24 | 0.296 | 2.156 | 1.5358 | 0.2504 |  P63/mcm
Si44 | 65 | 89 | 35 | 0.271 | 1.859 | 1.9308 | 0.0151 | Cmcm
Si${}_{\text{diamond}}$ | 91 | 151 | 62 | 0.223 | 1.472 | 2.3291 | 0.0000 | Fd-3m

The bulk modulus reflects the volumetric elasticity of the materials, while
Young modulus characterizes their tensile elasticity. Shear modulus, on the
other hand, signifies the material’s susceptibility to shearing when subjected
to opposing forces. The ratio of bulk modulus to shear modulus serves as an
indicator of the material’s brittleness or ductility. Specifically, if this
ratio falls below 1.75, the material is categorized as brittle, whereas values
exceeding this threshold denote ductility [43]. Furthermore, Poisson ratio can
also provide insights into the material’s brittleness or ductility, with
materials exhibiting a Poisson ratio exceeding 0.26 typically classified as
ductile [44].

Based on the information presented in Table 4, an analysis of the $B$/$G$ and
$n$ values for structures Si36(I), Si40, Si12(I), Si16(III), Si22, Si36(III),
Si16(IV), and Si44 suggests that these structures exhibit ductile properties.

Furthermore, when examining the $B$, $E$, and $G$ values, it becomes evident
that structures Si36(I) and Si40 have lower values compared to the other
structures. This can be attributed to the porous nature of these two crystal
structures.

Moreover, the relative enthalpy values presented in Table 4, which are
computed using the formula $\Delta
H=H/N-H_{\text{diamond}}/N_{\text{diamond}}$, where $H$ and $N$ represent the
enthalpy and the number of atoms in the respective structure, and
$H_{\text{diamond}}$ and $N_{\text{diamond}}$ correspond to the enthalpy and
the number of atoms in the silicon diamond structure, offer valuable insights
into the relative stability of these crystal structures. It is noteworthy that
Si12(III) emerges as the most stable structure, with Si20(II) and Si28 closely
following in terms of stability.

To ensure the stability of the newly developed structures, we undertook an
investigation focused on analyzing phonon scattering under varying pressures
(-5, 0.0, and 5 GPa). Initially, we allowed the structures to relax at two
distinct pressures: 5 and -5 GPa. Subsequently, we examined the phonon
dispersion of the structures for each pressure condition. The outcomes
revealed that all structures demonstrated stability under the specified
pressures, with the exception of Si16(III), which exhibited instability at 5
GPa pressure.

Figure 3: Energy of the structures as a function of volume.

Figure 3 depicts the energy variations of the structures in relation to volume
changes. As expected, the diamond silicon structure demonstrates the most
stability when compared to the other structures. Subsequently, Si12(III)
demonstrates a higher degree of stability in comparison to the remaining
configurations, securing the second-highest stability ranking. Following
closely behind is Si20(II), which attains the next most favorable stability
ranking.

In contrast, Si36(I) emerges as the least stable configuration, as indicated
by the graph’s data. Interestingly, Si40 exhibits relatively minimal energy
fluctuations in response to volume changes. According to Figure 3, it is
evident that structures Si40 and Si36(I) exhibit distinct characteristics
compared to other structures. These distinctions primarily arise from their
considerable volume, which is indicative of their porous nature. In order to
ascertain the stability of these two structures, a thorough investigation of
their phonon frequencies across various volumes was conducted. To achieve
this, we systematically generated eight distinct structures for each of the
aforementioned Si40 and Si36(I) structures, each featuring different volumes
within the range of 95% to 114% of the original structures’ volume.
Subsequently, we relaxed these structures to determine the equilibrium
positions of their constituent atoms under varying volumes and subsequently
computed their respective phonon frequencies. A notable discovery arising from
this investigation was the absence of any imaginary frequencies within the
calculated phonon spectra for these structures across the entire range of
volumes examined. This compelling result strongly supports the stability of
both the Si40 and Si36(I) structures.

###  3.3 Candidate Structures

Si12(III) stands out as the top choice for a solar cell when considering both
stability and GW band gap. However, if we focus solely on the GW band gap,
Si12(III), Si18(I), Si18(II), Si22, Si36(I), and Si40 are potential
candidates. Among these, Si40 appears promising for three reasons: I) Various
methods indicate a band gap of $\sim 1.35$ eV for Si40. II) The structure is
significantly different from Si diamond in the Born-Oppenheimer energy
landscape, making the transformation to Si diamond challenging due to a large
barrier. III) Si40 is ductile, making it practical for flexible solar cells.

###  3.4 Porous Structures with Void Tunnels

Figure 4: The ball-and-stick diagrams of the crystal structures.

In our research, we discover specific porous structures, like Si40,
characterized by void tunnels, setting them apart from others. These
structures exhibit hexagonal symmetry. To explore similar structures within
this symmetry, we employ an evolutionary algorithm. Our investigation leads us
to identify two stable structures (Si60 and Si24). Additionally, we attempted
to create Si52 by expanding upon the Si40 tunnel. These structures are
illustrated in Figure 4. Analyzing their properties, we found the HSE band
gaps for Si24, Si52, and Si60 to be 1.57, 1.00, and 0.76 eV, respectively.

These outcomes demonstrate a pattern: as the tunnel radius increases, the band
gap decreases. This phenomenon aligns with the concept of quantum confinement.
Thus, manipulating the tunnel radius offers a means to engineer and adjust the
band gap.

###  3.5 Machine Learning Calculations

The initial step in employing ML for predicting material properties involves
the preparation of input data. The generation of a suitable input derived from
a material has been a focal point of discussion among researchers. In general,
there are two basic approaches to construct the feature vector (or descriptor)
of a material. The first approach uses compositional information, typically
derived from the elemental and chemical properties and the stoichiometry of
the structure[45, 46, 47, 48]. The second approach relies on structural
information, such as atomic positions, exemplified by the Coulomb matrix[49]
and MBTR descriptors. However, since in this study, all the data contain the
same element (Si), the first approach is not applicable. This limitation
provides an excellent opportunity to assess the effectiveness of the second
approach. Additionally, we have generated descriptor for each structure from a
set of global structural properties, such as volume and energy of the
structures. This was accomplished by extracting 24 features from the outputs
of DFT calculations.

Furthermore, structural descriptors like MBTR and SOAP can be broadly
categorized into two main types: global and local descriptors [50, 38]. When
predicting a property of a material, it’s crucial to consider whether the
property is global or local in nature. This is important because attempting to
predict a global property using only local descriptors can lack physical
meaning. While it might be feasible mathematically, it may not make sense from
a physical perspective. If we intend to use a local descriptor to predict a
global property, such as band gap, we need to find a way to transform the
descriptor into a global one. This transformation is often achieved by
performing operations like inner or outer products on the feature vectors of
different local neighborhoods within the structure. As it detailed in the
literatures, the MBTR serves as a global descriptor, while the SOAP is a local
descriptor[50, 38]. Consequently, in this study, we have converted the SOAP
representation into a global representation to make it suitable for predicting
global properties like band gap.

In addition to these geometric descriptors, in this work, we use other
features obtained from the results of DFT calculations. The selected features
include the total energy (E), the total all-electron (AE) energy, internal
energy, one electron contribution, XC contribution, Ewald contribution,
Hartree (H) contribution,one-center PAW contribution, PAW Hartree energy of
AE, PAW Hartree energy pseudopotential (PS), PAW exchange-correlation (xc)
energy of AE, PAW xc energy of PS, total E-H with PAW, total E-XC with PAW,
Fermi energy, scf gap, the energy of the first three conduction bands at gamma
point denoted as ($\Gamma\\_E^{\prime}_{i}$) and the energy of the first four
valance bands at gamma point denoted as ($\Gamma\\_E_{i}$), along with volume
of the system.

Table 5 displays the Root Mean Squared Error (RMSE) values for band gap
predictions on both test and training data.

Table 5: The table presents Root Mean Squared Error (RMSE) and R2 values for predictions on both test and training data. These models were trained using data derived from SOAP and MBTR, and features extracted from the output of DFT calculations(DFT output). descriptor | ML algorithm | RMSE(train)eV | RMSE(test)eV |  R2(train) |  R2(test)
---|---|---|---|---|---
SOAP | KRR | 0.154 | 0.397 | 0.91 | 0.49
SVR | 0.039 | 0.353 | 0.91 | 0.49
NN | 0.183 | 0.393 | 0.87 | 0.45
CNN | 0.341 | 0.338 | 0.56 | 0.49
DT | 0.167 | 0.593 | 0.90 | 0.01
RF | 0.171 | 0.415 | 0.89 | 0.23
XGBoost | 0.113 | 0.403 | 0.95 | 0.27
MBTR | KRR | 0.105 | 0.403 | 0.96 | 0.25
SVR | 0.303 | 0.392 | 0.65 | 0.44
NN | 0.207 | 0.422 | 0.84 | 0.29
CNN | 0.285 | 0.299 | 0.69 | 0.60
DT | 0.099 | 0.623 | 0.96 | 0.01
RF | 0.136 | 0.343 | 0.93 | 0.59
XGBoost | 0.009 | 0.338 | 0.99 | 0.49
DFT output | KRR | 0.144 | 0.171 | 0.93 | 0.91
SVR | 0.105 | 0.169 | 0.96 | 0.90
NN | 0.050 | 0.215 | 0.99 | 0.81
CNN | 0.091 | 0.168 | 0.97 | 0.90
DT | 0.045 | 0.195 | 0.99 | 0.88
RF | 0.068 | 0.164 | 0.98 | 0.93
XGBoost | 0.011 | 0.155 | 0.99 | 0.91

As evident from the table, models trained on the features derived from DFT
calculations exhibit lower RMSE in band gap predictions. In the context of
this notation, "DFT output" signifies the features obtained through DFT
calculations. The R2(coefficient of determination) parameter values in the
table also indicate the superior performance of these models compared to those
trained on data obtained from MBTR and SOAP. The optimal models derived from
SOAP include the KRR and SVR, with RMSE values of 0.397 eV and 0.353 eV,
respectively, for the test data. Both models demonstrate R2 parameter values
of 0.91 for training data and 0.49 for test data. In the case of MBTR, the RF
and XGBoost models emerge as the best performers, with RMSE of 0.343 eV and
0.348 eV, respectively, for the test data. The R2 parameter values for these
models are 0.59 and 0.49, respectively, for the test data. When utilizing DFT
output features, the model accuracy significantly improves. The RMSE for all
models drops below 0.21 eV for test data, and the R2 parameter for test data
exceeds 0.80 across all models.

Figure 5: Band gap regression diagram using DFT output features and NN, SVR,
CNN, KRR RF, DT, and XGBoost ML models.

Figure 5 illustrates the relationship between band gap values predicted by ML
models (y-axis) and those obtained through DFTB calculations (x-axis). Only
models utilizing DFT output features are presented in this figure, while
graphs pertaining to models obtained from MBTR and SOAP are included in the
SI.

As evident from Figure 5, the models denoted as "DFT output + XGboost", "DFT
output + RF", and "DFT output + CNN" have emerged as the most effective
predictors.

Recognizing the significance of individual features within the context of ML
models enhances our ability to interpret the model’s performance and its
analytical outcomes. In this regard, the assessment of models utilizing
features extracted from DFT outputs, through the computation of Shapely
values[51] and the subsequent generation of beeswarm plots, serves as a
valuable analytical approach. Each model is accompanied by a beeswarm plot,
which effectively communicates the relative importance of each feature in the
prediction process. The graphical plots of Shapley values are available in SI.

The examination of these graphical representations unequivocally highlights
the feature labeled as "scf_gap" as possessing a superior degree of importance
when compared to the other features under consideration. This insight
underscores the pivotal role played by "scf_gap" in influencing the model’s
predictive performance.

Upon considering the second most crucial feature for band gap prediction, no
consistent ranking was observed across the various models. Specifically, in
the KRR and RF models, the "one_center_paw" feature claimed the second
position, whereas in the XGBoost and DT models, it occupied the third
position. Interestingly, in the NN, SVR, and CNN models, "one_center_paw" does
not appear among the top five important features.

A closer examination of feature importance from the third rank onward reveals
a distinct pattern in each model. There appears to be no predetermined order,
suggesting that different ML models do not uniformly prioritize features in
terms of their influence.

##  4 Conclusion

In this research, our primary objective was to predict and evaluate
crystalline silicon semiconductors, with a particular focus on their band gap
energy properties, with the aim of determining their suitability for
integration with solar cell technologies, thereby increasing the efficiency of
solar energy conversion systems. After extensive structural searches, we
screened the structure according to the band gap using fast DFTB calculations.
We recalculate band gaps using HSE and GW methods. The findings show that the
DFTB method with si-band SK is effective for screening band gaps and can be
applied in future studies. To ensure the uniqueness of these newly discovered
materials, a comprehensive analysis of their geometric configurations was
carefully performed. In addition, we performed a comprehensive evaluation of
the mechanical and dynamic stability of these structures. Our extensive
analyzes yielded compelling evidence of the dynamic stability of these silicon
structures, reinforcing their viability for integration into solar cell
technologies. We choose promising candidates, such as Si40 and Si12, with
potential applications in solar cells. Specifically, structures like Si40 show
void tunnels, suggesting not only an optimal band gap but also structural
flexibility. At the end, We used geometric properties and DFT-derived
properties to model and predict band gap energies. Subsequently, we performed
a critical assessment of the importance of selected features obtained from the
DFT outputs, and identified those elements that have the greatest impact on
the band gap prediction. In conclusion, our study provides a comprehensive
exploration of new direct and quasi-direct band gap silicon-based
semiconductor structures and positions them as promising candidates for next-
generation solar cell technology.

## References

  * Sreedevi  _et al._ [2019] P. Sreedevi, R. Vidya, and P. Ravindran, Solar Energy 190, 350 (2019).
  * Rühle [2016] S. Rühle, Solar Energy 130, 139 (2016).
  * Fan  _et al._ [2016] Q. Fan, C. Chai, Q. Wei, and Y. Yang, Phys. Chem. Chem. Phys. 18, 12905 (2016).
  * Zhang  _et al._ [2020] W. Zhang, C. Chai, Q. Fan, Y. Song, and Y. Yang, Results in Physics 18, 103271 (2020).
  * Zhao  _et al._ [2019] C.-X. Zhao, Y. Huang, J.-Q. Wang, C.-Y. Niu, and Y. Jia, Physics Letters A 383, 125903 (2019).
  * Wang  _et al._ [2014a] Q. Wang, B. Xu, J. Sun, H. Liu, Z. Zhao, D. Yu, C. Fan, and J. He, J. Am. Chem. Soc. 136, 125903 (2014a).
  * Xiang  _et al._ [2013] H. J. Xiang, B. Huang, E. Kan, S.-H. Wei, and X. G. Gong, Phys. Rev. Lett. 110, 118702 (2013).
  * Wei  _et al._ [2019] Q. Wei, W. Tong, B. Wei, M. Zhang, and X. Peng, Phys. Chem. Chem. Phys. 21, 19963 (2019).
  * Amsler  _et al._ [2015] M. Amsler, S. Botti, M. A. L. Marques, T. J. Lenosky, and S. Goedecker, Phys. Rev. B 92, 014101 (2015).
  * Botti  _et al._ [2012] S. Botti, J. A. Flores-Livas, M. Amsler, S. Goedecker, and M. A. L. Marques, Phys. Rev. B 86, 121204 (2012).
  * Lonie and Zurek [2011] D. C. Lonie and E. Zurek, Computer Physics Communications 182, 372 (2011).
  * Bahmann and Kortus [2013] S. Bahmann and J. Kortus, Computer Physics Communications 184, 1618 (2013).
  * Heydariyan  _et al._ [2018] S. Heydariyan, M. R. Nouri, M. Alaei, Z. Allahyari, and T. A. Niehaus, The Journal of Chemical Physics 149, 074313 (2018).
  * Seifert and Joswig [2012] G. Seifert and J.-O. Joswig, WIREs Computational Molecular Science 2, 456 (2012).
  * Elstner  _et al._ [1998] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, Phys. Rev. B 58, 7260 (1998).
  * Bartók  _et al._ [2010] A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Phys. Rev. Lett. 104, 136403 (2010).
  * Glass  _et al._ [2006] C. W. Glass, A. R. Oganov, and N. Hansen, Computer Physics Communications 175, 713 (2006).
  * Oganov  _et al._ [2006] A. R. Oganov, C. W. Glass, and S. Ono, Earth and Planetary Science Letters 241, 95 (2006).
  * Oganov  _et al._ [2010] A. R. Oganov, Y. Ma, A. O. Lyakhov, M. Valle, and C. Gatti, Reviews in Mineralogy and Geochemistry 71, 271 (2010), https://pubs.geoscienceworld.org/msa/rimg/article-pdf/71/1/271/2949866/271_Oganov_etal.pdf .
  * Sieck  _et al._ [2003] A. Sieck, T. Frauenheim, and K. A. Jackson, physica status solidi (b) 240, 537 (2003).
  * Aradi  _et al._ [2007] B. Aradi, B. Hourahine, and T. Frauenheim, The Journal of Physical Chemistry A 111, 5678 (2007), pMID: 17567110, https://doi.org/10.1021/jp070186p .
  * Bartók  _et al._ [2018] A. P. Bartók, J. Kermode, N. Bernstein, and G. Csányi, Phys. Rev. X 8, 041048 (2018).
  * Bartók  _et al._ [2013] A. P. Bartók, R. Kondor, and G. Csányi, Phys. Rev. B 87, 184115 (2013).
  * Bartók-Pártay [2010] A. Bartók-Pártay,  , 5 (2010).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Hedin [1965] L. Hedin, Phys. Rev. 139, A796 (1965).
  * Becke [1993] A. D. Becke, The Journal of Chemical Physics 98, 1372 (1993), https://pubs.aip.org/aip/jcp/article-pdf/98/2/1372/11046762/1372_1_online.pdf .
  * Markov  _et al._ [2015] S. Markov, B. Aradi, C.-Y. Yam, H. Xie, T. Frauenheim, and G. Chen, IEEE Transactions on Electron Devices 62, 696 (2015).
  * Marini  _et al._ [2009] A. Marini, C. Hogan, M. Grüning, and D. Varsano, Computer Physics Communications 180, 1392 (2009).
  * Mostofi  _et al._ [2014] A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer Physics Communications 185, 2309 (2014).
  * Corso [2023] A. D. Corso, “Thermo pw,”  (2023).
  * Prandini  _et al._ [2018] G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet, and N. Marzari, npj Computational Materials 4, 72 (2018), http://materialscloud.org/sssp.
  * Eriksson  _et al._ [2019] F. Eriksson, E. Fransson, and P. Erhart, Advanced Theory and Simulations 2, 1800184 (2019).
  * Abadi  _et al._ [2015] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,”  (2015), software available from tensorflow.org.
  * [35] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, Journal of Machine Learning Research 12, 2825.
  * Huo and Rupp [2022] H. Huo and M. Rupp, Machine Learning: Science and Technology 3, 045017 (2022).
  * Marc O. J. Jäger [2018] F. F. C. L. H. A. S. F. Marc O. J. Jäger, Eiaki V. Morooka, npj Computational Materials 4, 2057 (2018).
  * Himanen  _et al._ [2020] L. Himanen, M. O. Jäger, E. V. Morooka, F. Federici Canova, Y. S. Ranawat, D. Z. Gao, P. Rinke, and A. S. Foster, Computer Physics Communications 247, 106949 (2020).
  * Valle and Oganov [2010] M. Valle and A. R. Oganov, Acta Crystallographica Section A: Foundations of Crystallography 66, 507 (2010).
  * Togo and Tanaka [2018] A. Togo and I. Tanaka, arXiv preprint arXiv:1808.01590 5 (2018), https://arxiv.org/abs/1808.01590.
  * Wang  _et al._ [2014b] Q. Wang, B. Xu, J. Sun, H. Liu, Z. Zhao, D. Yu, C. Fan, and J. He, Journal of the American Chemical Society 136, 9826 (2014b), pMID: 24971657, https://doi.org/10.1021/ja5035792 .
  * Lee  _et al._ [2014] I.-H. Lee, J. Lee, Y. J. Oh, S. Kim, and K. J. Chang, Phys. Rev. B 90, 115209 (2014).
  * Pugh [1954] S. Pugh, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 45, 823 (1954), https://doi.org/10.1080/14786440808520496 .
  * J. J. Lewandowski * and Greer [2005] W. H. W. J. J. Lewandowski * and A. L. Greer, Philosophical Magazine Letters 85, 77 (2005), https://doi.org/10.1080/09500830500080474 .
  * Zhang and Ling [2018] Y. Zhang and C. Ling, npj Computational Materials 4 (2018), 10.1038/s41524-018-0081-z.
  * Legrain  _et al._ [2017] F. Legrain, J. Carrete, A. van Roekeghem, S. Curtarolo, and N. Mingo, Chem. Mater. 29 (2017), 10.1021/acs.chemmater.7b00789.
  * Jha  _et al._ [2018] D. Jha, L. Ward, A. Paul, W.-k. Liao, A. Choudhary, C. Wolverton, and A. Agrawal, Scientific Reports 8 (2018), 10.1038/s41598-018-35934-y.
  * Jain and Bligaard [2018] A. Jain and T. Bligaard, Phys. Rev. B 98, 214112 (2018).
  * Rupp  _et al._ [2012] M. Rupp, A. Tkatchenko, K.-R. Müller, and O. A. von Lilienfeld, Phys. Rev. Lett. 108, 058301 (2012).
  * Damewood  _et al._ [2023] J. Damewood, J. Karaguesian, J. R. Lunger, A. R. Tan, M. Xie, J. Peng, and R. Gómez-Bombarelli, Annual Review of Materials Research 53, 399 (2023), https://doi.org/10.1146/annurev-matsci-080921-085947 .
  * Lundberg and Lee [2017] S. M. Lundberg and S.-I. Lee, in  _Advances in Neural Information Processing Systems_, Vol. 30, edited by I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Curran Associates, Inc., 2017).