Materials discovery via CALYPSO methodology

It is well known that for the ordered materials the atomic packings are governed by the symmetry rules (e.g. determined by space groups for three-dimensional (3D) solids or point groups for zero-dimensional (0D) isolated clusters or molecules). For structure prediction, the implementation of symmetry constraints is crucially important since it can reduce the structure variables to be optimized, and thus significantly reduce the searching space of PES. Taking a simple example of crystals with a monoclinic symmetry, compared with the case excluding the symmetry constraints, the number of variables of lattice parameters to be optimized is reduced from 6 to 4 (there are two fixed lattice angles of 90°), and the varying range of fractional atomic coordinates is reduced to 0–0.5 other than 0–1.0 due to the presence of mirror operation.

In our CALYPSO method, different strategies dealing with 3D crystals, two-dimensional (2D) layers or surfaces, and 0D-isolated systems are designed to randomly generate candidate structures. Particularly, the generation of random structures is constrained by 230 space groups for bulk crystals, the 17 in-plane space groups for 2D layers or surfaces, and the 32 point groups for isolated systems (e.g. molecules, clusters, and nanoparticles). Random generation of candidate structures with symmetry constraints in structure prediction was pioneered by several groups [9, 16, 17], and it was later adopted in [25].

During structure searches, a large number of trial structures residing in the same basin of PES will be generated, each of which is then subjected to local structure optimization. Repetitive optimizations of these structures will significantly slow down the search efficiency. As such, it is highly demanding to have a technique that can unbiasedly fingerprint a structure and measure the similarity between involved structures. With the help of such technique, one can deliberately remove duplicated structures to avoid the waste of computational time on their structure optimization. By taking advantage of structural bond information (bond lengths and angles), we have developed a structural characterization technique named as bond characterization matrix (BCM) and implemented it into our CALYPSO method. This is realized by constructing a set of modified bond-orientational order metrics [26] for structure quantification, where spherical harmonic and exponential functions were used to characterize the bond angles and bond lengths, respectively. Specifically, for a given structure, a bond vector $\overrightarrow {r_{ij}}$ between atoms i and j was defined if the interatomic distance is less than a given cutoff distance. $\overrightarrow {r_{ij}}$ is associated with the spherical harmonics Y_lm(θ_ij, φ_ij), where θ_ij and φ_ij are the polar angles, and the weighted average over all bonds formed by for instance A and B atoms can be derived by the equation:

$$\begin{equation} \bar{{Q}}_{lm}^{\delta_{AB}} =\frac{1}{N_{\delta_{AB}}}\sum\limits_{i\in A,j\in B} {{\rm e}^{-\alpha (r_{ij} -b_{AB} )}Y_{lm} (\theta_{ij} ,\phi_{ij} )}, \end{equation} \tag{ 1 }$$

where δ_AB and $N_{\delta_{AB}}$ denotes the type and the number of bonds, respectively. Only even-l spherical harmonics are used in equation (1) to guarantee the invariant bond information with respect to the direction of the bonds. In order to avoid the dependence on the choice of reference frame, it is important to consider the rotationally invariant combinations [26],

$$\begin{equation} Q_{l}^{\delta_{AB}}=\sqrt {\frac{4\pi}{2l+1}\sum\limits_{m=-l}^l {\left| {\overline Q_{lm}^{\delta_{AB}}} \right|^{2}}} , \end{equation} \tag{ 2 }$$

where each series of $Q_{l}^{\delta_{AB}}$ for l = 0, 2, 4, 6, 8, 10 can be used to represent a type of bonds, thus being an element of BCM. By including the complete structure information on bond lengths and bond angles, this technique has the capability of unambiguously fingerprinting structures. As a result, the similarity of two structures can be quantitatively evaluated by the Euclidean distance between their BCMs,

$$\begin{equation} D_{uv} =\sqrt {\frac{1}{N_{type}}\sum\limits_{\delta_{AB}} {\sum\limits_l {(Q_{l}^{\delta_{AB} ,u} -Q_{l}^{\delta_{AB} ,v} )^{2}}}} , \end{equation} \tag{ 3 }$$

where u and v denote two individual structures and N_type is the number of bond types. Currently, the BCM technique has been applied to the structure prediction for 3D crystals [24] and 0D clusters [27].

Swarm intelligence [28] algorithm inspirited by natural biological systems (e.g. ants, bees or birds) is based on the collective self-organized behaviors of individual agents following simple and 'intelligent' rules. Typical swarm intelligence schemes include PSO [29], ant colony system [30], artificial bee colony [31], and so on [32, 33]. The swarm intelligence algorithm has been applied to a variety of fields in engineering and chemical science [34].

As a stochastic global optimization algorithm, PSO was proposed by Eberhart and Kennedy in the mid 1990s [37], inspired by the social behavior of birds flocking or fish schooling. Within the framework of PSO, the behavior of each individual is affected by both the local best individual and the global best individual. Moreover, an individual can learn from its past experiences to adjust its flying speed and direction. Therefore, all the individuals in the swarm can quickly converge to the global position. PSO was attempted for isolated systems (small clusters and molecules) by Call et al in 2007 [35]. We independently initialized the idea on application of PSO into structure prediction and introduced the PSO algorithm for the first time into extended 2D and 3D systems. Within our CALYPSO method, structures are evolved in the PES through the defined velocity of individual following the formula (4). The new velocity of each structure v^t+1 is calculated based on its previous velocity (v^t), previous location (x^t), current location (pbest^t) corresponding to an achieved best fitness value of this individual, and the global location (gbest^t) corresponding to the best fitness for the entire population according to equation (5).

$$\begin{equation} x^{t+1}=x^{t}+v^{t+1}, \end{equation} \tag{ 4 }$$

$$\begin{equation} v^{t+1}=\omega v^{t}+c_{1} r_{1} ({\rm pbest}^{t}-x^{t})+c_{2} r_{2} ({\rm gbest}^{t}-x^{t}). \end{equation} \tag{ 5 }$$

Here ω denotes the inertia weight, which is dynamically varied and decreases linearly from 0.9 to 0.4 during the iteration; c₁ and c₂ are the self-confidence factor and swarm confidence factor, respectively; r₁ and r₂ are random numbers, which are distributed in the range of [0, 1]. The velocity updating (equation (5)) includes two random parameters r₁ and r₂ to ensure good coverage of the searching space and avoid undesired entrapment in some local minimum. By making use of the swarm intelligence and by self-improving structures, PSO is best known for its ability to conquer large barriers of free energy landscapes.

Currently, two versions of PSO have been implemented in CALYPSO method: the global version and the local one. The schematic diagram of global PSO is illustrated in figure 1(a). In this version, there is only one global best structure acting as the attractor for the entire structure population and all particles seek new positions only in the regions close to the uniquely overall best position. The global PSO has the advantage of being able to reach convergence fast for small systems (i.e. the number of atoms in the simulation cell smaller than 30). But it may be less effective for the larger systems containing more than 30 atoms per simulation cell, where the free energy landscape becomes much more complex [36]. The problem can be attributed to the loss of diversity owing to the existence of only one attractor for the entire population. This deficiency can be remedied by introducing the local PSO scheme, in which the information is diffused into small parts of the swarm [37].

Zoom In Zoom Out Reset image size

The schematic diagram of the local PSO is illustrated in figure 1(b). In this version, each particle selects a set of other particles as its neighbors and its velocity is adjusted according to both its position and the best position achieved so far in the community formed by its neighborhood. Thus at each iteration, the particle will move towards its own best position and the best position of its local neighborhood, instead of the overall best position in the swarm. By maintaining multiple attractors that lead to the formation of different structural motifs, the local PSO allows an unbiased and efficient exploration of larger space on PES, thus effectively avoiding the premature convergence during structure searches.

Maintaining high structural diversity is critically important for ensuring the success of population-based evolutionary structure searching method. When the loss of structural diversity occurs during the structural evolution, the search is likely to become stagnant at one predominated local minimum, thus leading to the premature convergence [38]. To avoid this problem, we have included a certain percentage of new random structures per generation during structural evolution in an effort to enhance the structural diversity. In this technique, the symmetry is used to differentiate structures and the appearance of the structures with identical symmetry is forbidden at a certain (∼80%) probability. This allows the samplings to cover different regions of the search space and thus forces the search to explore new minima. It has been proved that this technique can significantly enrich the structural diversity, thus effectively reducing the stagnation rate of CALYPSO method [24].

The PES of a material can be regarded as a multi-dimensional system having a large number of hills and valleys with saddle points connecting them. Each of the valleys forms a local basin of attraction on the PES. To get to the local minimum of basins (figure 1(c)), we can use structural optimization techniques by minimizing the total energy of systems (e.g. steepest descents [39], conjugate gradient algorithm [40] or Broyden–Fletcher–Goldfarb–Shanno algorithm [41]), after which the free energy of candidate structures is evaluated, and used as the fitness function of evolution. Local structural optimization increases the computational cost on each individual, but effectively reduces the noise of energy landscape and produces physically reasonable structures for further evolution. As such, local structural optimization acts as a crucial step for all the modern structure search methods. Our CALYPSO method currently has interfaces with a variety of ab initio and force-field based total energy packages (e.g. VASP [42], SIESTA [43], Quantum-Espresso [44], CASTEP [45], CP2K [46], Gaussian [47], DFTB+ [48], LAMMPS [49] and GULP [50]) for local structural optimization. Other external total energy programs can also be easily interfaced with CALYPSO at users' request.

In crystallography, a crystal structure can be regarded as an infinitely repeating array of 3D boxes, known as unit cells (i.e. Bravais lattice + basis). There are two types of variables to define the unit cell: lattice parameters (the lengths of three lattice vectors and three angles among them) and atomic fractional coordinates (three components coded as the fractions of lattice vectors for each atom). For structure searches of 3D crystals, both of two types of variables are evolved simultaneously to achieve the global minimum on the PES.

Most of known structure prediction methods/codes have been successfully applied to the investigation of 3D solids. The detailed descriptions of these methods and some of their applications have been presented in our recent perspective article [23]. It is not feasible to give a precise ranking on searching efficiency among these different codes since they are all stochastic methods and their performance is sometimes dependent of systems. Fortunately, different methods can give the same converged results for most of calculations as exemplified by the predicted identical high pressure structures of N₂ [51], Li [52, 53] and Xe–Fe systems [54] as derived by CALYPSO [17] and AIRSS [16] methods/codes.

In order to deal with the 2D layered materials, which consist of planar (or buckled) monolayer or multiple layers, we implemented a 2D structure search module in CALYPSO method. As illustrated in figures 2(a) and (b), the configuration space of 2D layered structures contains two regions: the layered material part and the vacuum part. The latter is used to ensure that the studied layer materials are isolated from its nearest-neighboring periodic image. In addition to these, a distortion parameter perpendicular to the in-plane layer (Δz as denoted in figure 2(a)) can be turned on to search for the buckled layer, and a van der Waals gap parameter (i.e. the distance between two adjacent layers as in figure 2(b)) is introduced for multiple layered systems.

Zoom In Zoom Out Reset image size

To perform structure searches for surfaces, a slab module is adopted, where the surface is simulated by a thick enough film separated by a vacuum along the direction perpendicular to the surface (figure 2(c)). The thin film consists of three regions: the bulk material region, the unreconstructed surface and the reconstructed surface. Since the surface properties usually show quick convergence along the perpendicular direction, for most crystalline surfaces 6–8 atomic layers are enough to represent the bulk region. While the bulk region keeps fixed to reserve the bulk nature of materials, the atoms in the unreconstructed surface region (usually composed of 2–4 layers) are subject to local structural relaxation, but they are not involved in the structural evolution. The atoms in the reconstructed surface region are fully evolved during structure search, and the choice of thickness for this region depends strongly on specific energetic situation of systems. The fitness function during structural evolution is set to the surface excess free energy (that represents the energetic stability of surface reconstructions) following the equation

$$\begin{equation} \Delta \gamma (\mu )=\frac{1}{A}\left( {E_{{\rm surf}}^{{\rm tot}} -E_{{\rm ideal}}^{{\rm tot}} -\sum\limits_i {n_{i} \ast \mu_{i}}} \right), \end{equation} \tag{ 6 }$$

where A is the area of studied surface (usually in the unit of 1 × 1 planar unit cell), $E_{{\rm surf}}^{{\rm tot}}$ and $E_{{\rm ideal}}^{{\rm tot}}$ are the total energies of reconstructed and ideal (unreconstructed) surfaces, n_i and μ_i are the number and chemical potential of the ith type specie in the surface region.

Clusters (or nanoparticles) belong to 0D non-periodic material systems, resembling the large molecules in chemistry. The CALYPSO method has also been generalized to carry out structure search for such non-periodic systems (see (figure 2(d)) for their structural models). Different from the periodic system with translational symmetry, in isolated systems only the point group symmetries were utilized in generating candidate structures. To comply with usual total energy calculations requiring periodic boundary condition, a big box is built, where the cluster is located at its center. The vacuum surrounding the clusters should be large enough to avoid the interaction between the cluster and its periodic images. For such non-periodic calculations, Cartesian coordinates are used straightforwardly to represent structures. The module of cluster structure search has been extensively benchmarked in the Lennard–Jones cluster systems with varied sizes (up to 150 atoms) [27].

Stoichiometry is one of the fundamental thermodynamic variables for determining materials' properties. It is also covered in CALYPSO code by implementing a module of automatic variation of chemical compositions for structure search. With the aim of identifying at one time all the energetically stable structures throughout the whole composition range, this module has the capability of constructing a full convex hull of formation energy versus composition for a given system. The schematic of the convex hull evolution during structure search is illustrated in figure 2(e): at the nth step, the available lowest–energy structures (i.e. a, b, c and d) at different stoichiometries are used to construct an approximate convex hull. Afterward, the new convex hull is updated by considering the newly emerged structure a', leading to the formation of a new hull through a', c and d at (n + 1)th step. The structure b is reasonably discarded since its formation energy lies above the tie line between the neighboring pair of structures a' and c. The direct search for convex hull by using this module is usually more efficient than the study involving separate structure searches at different compositions.

In some cases, partial structural information (such as unit cell, partial atomic positions, molecular unit or structural motif) is available from experiments. Taking full advantage of such information is extremely useful for structure search because it can significantly reduce the configuration space of candidate structures. For this we have implemented a module flexibly using known partial structural information into CALYPSO method. For example, in the system with the experimentally known molecular unit, the Z-matrix method [55] is employed to represent the rigid molecular unit, which can only move or rotate as a whole unit during structural evolution. The prior knowledge of other structural information such as lattice parameters, partial atomic positions or space group can also be used to accelerate structure search and guarantee the search proceeding along the right direction.

The conventional way of material research usually starts from determining the structure of a material just synthesized, then followed by systemically characterizing its properties. The shortcoming of this approach is the absence of predictable power, i.e. one can never know what material structure can offer some desired property. One promising way to overcome this shortcoming is the computational design of functional materials with desired target properties by using advanced supercomputing facilities. This can be readily done by replacing the fitness function of total energy with various material properties. In this respect, we developed a material design module based on CALYPSO method to design new functional materials. For instance, by choosing hardness or band gap as the fitness function, CALYPSO can be used to design superhard [56] and optoelectronic materials [57, 58], respectively. Of course, for such functionality-driven design of materials, the criterion of energetic stability has to be properly considered.

Electrides, in which localized electrons occupy interstitial regions of the crystal and behave as anions, while the ionic cores are cations, are known in various organic compounds [59]. Until 2009, Ma and coworkers has first introduced this concept to describe the transparent insulating phase of sodium at high pressure >200 GPa [60]. Since then, a series of high-pressure elemental electrides [52, 61–64] have been theoretically proposed. There, the CALYPSO method plays an important role in discovery of elemental electrides such as Mg [64] and Li [52].

3.1.1.1. Magnesium.

3.1.1.2. Lithium.

Lithium exhibits complex structural behavior at high pressure, which leads to topics of interest in its structural transitions. At ambient condition, it adopts a bcc structure and transforms to a complex cubic I-43d structure at 42 GPa via two intermediate phases (fcc and rhombohedral R-3 m) [68]. Above the stable pressure range of I-43d structure (>70 GPa), experimental measurements [69–72] suggest the existence of semiconducting broken-symmetry phase. However, its crystal structure remains mystery. Since then, much effort [62, 73–76] has been devoted to the determination of this insulating phase. In 2011, we extensively explored the energy landscape of dense Li using the CALYPSO method and predicted a novel insulating structure of Aba2-40 (figure 3(a)) to be most stable at 60–80 GPa [52]. Our calculations reveal that this insulating electronic state emerges in the Aba2-40 phase because of core exclusion and the localization of valence electrons in the voids of crystal, which results in the formation of a new high-pressure electride similar to that in transparent Na [60]. This theoretical prediction was confirmed by the independent high-pressure experiments [63] and another theoretical calculation [53].

Zoom In Zoom Out Reset image size

3.1.2.1. Hydrogen-rich superconducting compounds.

3.1.2.2. Hydrogen-free superconducting compounds.

High energy density material (HEDM) becomes increasingly important for the development of clean and efficient energy material. High-nitrogen compounds are promising candidates for HEDM because of their remarkable high positive formation heats. Among these compounds, polymeric nitrogen and alkali metal azides have attracted growing interest due to the high nitrogen content.

3.1.3.1. Polymeric nitrogen.

3.1.3.2. Alkali metal azides.

3.1.4.1. Light-elements compounds.

3.1.4.2. Transition-metal light-element compounds.

Li is the lightest metallic element in periodic table. Li clusters are thus considered to be prototype systems for understanding the various physical properties of simple metal clusters. Since the seminal experimental work by Knight et al [127], numerous theoretical studies [128–132] have been performed to understand the structures of Li clusters. Geometric structures for Li clusters with sizes up to 147 atoms have been proposed based on the density functional calculations. As the first realistic application of the CALYPSO method on cluster structure prediction, medium-sized Li_n(n = 20, 40, 58 [27], which are magic numbers according to the jellium model) clusters have been revisited. For Li₂₀ cluster, the predicted lowest energy structure is in agreement with the previous results [133]. For Li₄₀ and Li₅₈, new putative global minima with polyicosahedral and Mackay icosahedral motifs have been proposed [27].

The discovery of C₆₀ fullerene [134] had opened up a new field of fullerene research. Since then, the search for fullerene-like structures formed by elements other than carbon has become one of the most challenging tasks in physical, chemical, and material sciences. Boron, the early neighbor of carbon in the periodic table, is also able to form sp² and sp³ bonding and thus is believed to be the best candidate after carbon to form fullerene-like structure. Recently, CALYPSO structure searching of 38-atom boron cluster leads to the discovery of an all-boron B₃₈ fullerene analogue (figure 3(d)). The structure is unexpectedly chemically stable, as manifested by a large energy gap and a high double aromaticity. Moreover, by means of the CALYPSO structure searching and high temperature molecular dynamics simulations, the structures of 55-atom transition metal cluster were explored by Li et al [135]. They found that fcc- or hcp-fragment structures are more stable than the icosahedral structure and thus the widely perceived magic size of 55 is shifted to its nearby even numbers.

Recently, the intriguing planar wheel-type geometries of ${\rm B}_{8}^{2-}$ and ${\rm B}_{9}^{-}$ clusters have led to the discovery of a series of transition-metal-centered borometallic wheels: ${\rm FeB}_{8}^{-}$, ${\rm FeB}_{9}^{-}$, ${\rm CoB}_{8}^{-}$, ${\rm RuB}_{9}^{-}$, ${\rm RhB}_{9}^{-}$, ${\rm IrB}_{9}^{-}$, ${\rm NbB}_{10}^{-}$ and ${\rm TaB}_{10}^{-}$, where the ${\rm NbB}_{10}^{-}$ and ${\rm TaB}_{10}^{-}$ hold the highest coordination number in planar species [136]. Based on the electronic and geometric principles, Liao et al [137] have designed a series of boron wheels enclosing planar hypercoordinate M atoms (M is a 2nd or 3rd period element), then they performed CALYPSO structure searches to examine the stabilities of these newly designed clusters. It is found that most of these boron metallic wheels are only local, rather than global minima. However, the BeB₈ triplet planar wheel has been revealed to be a new member of this planar hypercoordinate family.

CALYPSO method for 2D structure prediction has been benchmarked by several 2D systems including elemental (C and B) or binary compounds (BN). All of the experimentally and theoretically known structures are successfully reproduced by CALYPSO method, validating our approach [138, 139]. Using CALYPSO method, we studied the structure of a new family of layered B_x N_y compounds (B₂N₃, B₃N₄ and B₃N₅) and uncovered a serial of 2D structures with different structural motifs. Particularly, an unexpected novel B₃N₅ 2D sheet (figure 4(a)) formed by entirely odd-numbered ring is firstly found for a planar compound [138]. Electronic structure calculation shows that it is a semiconductor with a band gap ∼2.2 eV, indicating its potential applications as electronic and optical materials. In addition, a number of novel 2D sheet structures [140–150] have also been designed by CALYPSO method.

Zoom In Zoom Out Reset image size

Mechanically assembled 2D multi-layer materials are expected to show combined functionality of the individual layers. Recently, the first-principles calculations [151] show that the band gap and effective electron mass of graphene/BN heterostructures might be tuned by the interlayer spacing and stacking arrangement of the individual layers. Thus, 2D multi-layer materials have potential applications for the development of novel technologies. Using CALYPSO method, we have designed the structures of the C/Si/BN sandwich material. Two tri-layer sandwich structures (S1 and S2) formed by stacking of graphene, trigonal Si and hexagonal BN (h-BN) layers were proposed. As shown in figures 4(b) and (c), they share the same stacking arrangement of C and Si layers. Specifically, two distinct C atoms locate above the Si atom and the center of trigonal ring formed by Si atoms, respectively. The difference between these two structures is the stacking arrangement of h-BN. For S1 structure, the B and N atoms locate at center of Si trigonal ring, while the B atom sits below the Si atom for S2 structure. The effects of stacking sequences and arrangements for design of novel multi-layer materials can be explored in the future.

There are numerous exciting opportunities in design of 2D materials with novel functionalities where CALYPSO can play an important role in such designs. These researches will inevitably lead to the discovery of novel materials with unanticipated properties.

Diamond (1 1 1) and (1 0 0) surface are of both fundamentally and technically interest. From a fundamental point of view, it is interesting to understand the driving forces of the reconstruction behavior and bonding feature, especially compared with other group-IV semiconductors, i.e. silicon and germanium surfaces [152–154]. On the other hand, the diamond (1 0 0) and (1 1 1) surface are the two dominant growth planes in the production of diamond via chemical vapor deposition method [155]. Using CALYPSO surface reconstruction search method [156], we successfully identified the previously reported diamond (1 1 1) Pandey-chain model (figure 5(a)) and (1 0 0) dimer reconstruction model (figure 5(b)). The diamond (1 1 1) surface exhibits a dramatic reconstruction to the Pandey-chain model. The upmost carbon atoms form zigzag π-bonded chains that run in [0 1 1] direction. Besides, as shown figure 5(c), a hitherto unexpected surface reconstruction of diamond (1 0 0) featuring self-assembly of carbon nanotubes (CNTs) arrays has been found by CALYPSO surface reconstruction prediction module [156]. The surface with self-assembled CNTs is energetically degenerate to the dimer structure at normal conditions, but become much favorable under a small compressive strain or at elevated temperatures [156].

Zoom In Zoom Out Reset image size

Hydrogen is generally necessary for diamond film growth via the chemical vapor deposition method. Furthermore, hydrogenated diamond (1 0 0) surface can deliver negative electron affinity [157, 158], and high p-type conductivity [159]. Using CALYPSO method, we successfully reproduced the hydrogenated dimer surface reconstruction. As shown in figure 5(d), each surface carbon atom bonds to a hydrogen atom, forming four-fold sp³ hybridization. We also found that diamond (1 0 0) surface will form 2 × 2 with trough surface reconstruction (figure 5(e)) under hydrogen-rich environments.

Recently, we have extended CALYPSO method into inverse design of superhard materials [56]. In this method, the hardness was adopted as the fitness function, other than the traditional structure prediction in seeking for global energy minimum. By combination with the first-principles energetic calculation, hardness versus energy map can be constructed for a given chemical system, from which the energetically favorable superhard structures are readily accessible [56]. Using this approach, we designed various superhard materials among three typical superhard light-elements systems including elemental carbon, binary B–N, and ternary B–C–N compounds [56]. Nearly all the experimentally known and earlier theoretical superhard structures have been successfully reproduced, indicating that our approach is reliable and can be widely applied into design of new superhard materials. Besides the known structures, our method predicted several new meta-stable superhard structures. For elemental carbon, two orthorhombic meta-stable structures with symmetry of P222₁ and Imma were found to be superhard materials, whose theoretical hardness values are comparable to that of diamond. These two phases might be synthesizable from compressed graphite since the transition pressures of graphite into them are lower than that into M-carbon [160], which has been experimentally synthesized [161].

Silicon plays an important role in the field of photovoltaic energy production. Most of first-generation solar cells use bulk silicon as the absorber layer [162]. However, the indirect gap nature of diamond Si with a large energy difference between the direct gap (3.4 eV) result in the low efficient absorption of sunlight [163]. Thus, it is desirable to discover new Si phases with suitable direct gap. Recently, CALYPSO was adopted to explore the new phases of silicon with direct bang gaps, where the fitness function is the band gap, rather than the total energy of the structure. A cubic Si phase with a quasi-direct gap of 1.55 eV was designed [57]. Later, more metastable allotropes of silicon with direct or quasi-direct band gaps of 0.39−1.25 eV were predicted by CALYPSO method [164]. These structures can absorb sunlight with different frequencies, providing appealing features for application in the tandem multi-junction photovoltaic modules.