Learning the stable and metastable phase diagram to accelerate the discovery of metastable phases of boron

Boron is widely used in a variety of applications, for example, in superconducting materials,¹ in glasses,² and in the semiconductor industry as a p-type doping element and as dielectrics.³ In addition, boron is also a common component in many superhard materials, such as WB₄⁴ and boron nitride.⁵ The wide applicability of boron and its tendency to form clustered compounds with varying coordination numbers have attracted the attention of many researchers to compute the phase diagram of boron.^6–9 Computational search for the equilibrium phase diagram of compounds focuses on the thermodynamic aspect and searching for structures with the lowest energy. However, many phases are often trapped in local minima in free energies or metastable states. Thus, an equilibrium phase diagram does not give an overall picture of the competing metastable phases, which provides a richer context on the phases formed and their synthesizability. To attain the target synthesizability of these metastable structures, computation of both the ground state phase diagram and the metastable phase diagram, as a function of pressure and temperature, becomes mandatory. In addition, many phases with exceptional properties, such as hardness and bandgap, could be metastable, such as superhard boron T-50 (T-B₅₀).⁴ Thus, computing and constructing the metastable phase diagram would give researchers an estimate of the synthesizability of the metastable phases.

The metastable phase diagram is a 2D slice of the complex 3D free energy surfaces of the various competing phases as a function of temperature and pressure. An initial challenge in the construction of the metastable phase diagram is to utilize an efficient algorithm for structure optimization that can identify both global (ground state) and local (metastable) minima in the energy landscape of the configuration space. The subsequent challenge involves mapping the free energy surface, or equation of state, for each metastable phase as a function of key thermodynamic variables (such as pressure, temperature, and composition) within the desired phase information range. However, this step can often become computationally demanding when dealing with numerous metastable configurations, necessitating the use of a surrogate model to approximate the costly calculations of a first-principles based method, such as ab initio molecular dynamics. Once the equation of state for all phases is computed, a final challenge is to classify and identify phase boundaries and domains of metastable equilibrium, referring to areas in the phase diagram where a metastable structure remains distinct from lower energy structures. A new procedure that overcomes these challenges for estimating the phase diagram and the metastable phase diagrams at different energies from ground state has recently been established by us and studied for a carbon system.¹⁰ In this procedure, the temperature–pressure dependent Gibbs free energies of the competing phases are estimated using a quasi-harmonic approach to estimate entropy. The stable and metastable phases under each temperature–pressure condition are calculated based on the Gibbs free energies of the phases under those conditions.

Computational search for different metastable phases of materials has always been a challenging task because of the high dimensionality of the search space. Evolutionary crystal structure prediction algorithms, such as USPEX,¹¹ CALYPSO,¹² and XTALOPT,¹³ have been implemented successfully to search for ground state structures. Machine learning methods have tried to bridge this gap by using a high throughput approach^14–17 to sample material properties cheaply. Recently, a reinforcement learning coupled with tree methods, Monte Carlo Tree Search (MCTS),^18–23 has emerged as a powerful algorithm to explore search spaces, such as in game spaces and drug discoveries.^24,25 MCTS is efficient due to its relatively low computational cost and its balance in exploring the space and exploitation of the already searched space, i.e., the balance in exploration–exploitation, and is explained in detail later. MCTS has been adapted for structure search and has found success in predicting optimal defect configurations in 2D materials²⁰ and phase transitions in phosphorene,²⁶ in predicting bond order potentials for clusters,²⁷ and recently, in Crystal Structure Prediction (CSP) for materials of varying classes and dimensionality.²⁸ Inspired by its success, we use MCTS to search the configurational–compositional space for boron structures and use the identified configurations to construct the equilibrium and metastable phase diagrams of boron. In this work, we report on some new metastable boron structures sampled by the MCTS approach. We construct the equilibrium phase diagram with the sampled boron structures and make a comparison with existing reports. Next, the metastable phase diagram is constructed to show the temperature dependence of select competing boron phases, which were sampled using MCTS. Slices of metastable phase diagrams at different free energy levels from ground state equilibrium are obtained, and the equations of state are compared to give an idea of the synthesizability of these phases.

Density functional theory calculations were performed using the Vienna Ab initio Simulation Package (VASP),^29,30 employing periodic boundary conditions, a plane wave basis set, the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA) exchange–correlation functional,³¹ and the projector-augmented wave method.³² All computational parameters were chosen such that the calculated total energy differences were converged to within 1 meV/atom. This included a 500 eV cutoff energy for the plane wave basis set expansion and a Γ-centered k-point grid with a k-point density of 2000.

The structure search is performed by implementing the Monte Carlo Tree Search (MCTS) methodology as implemented in the CASTING package.²⁸ The algorithm optimizes the structures in the configurational space (a, b, c, α, β, and γ) and compositional space (n_atoms). An upper limit of 50 atoms is set as the maximum number of atoms reached in the compositional space. The modules of the CASTING framework are initialized to optimize the objective function using the Upper Confidence of Bounds (UCB) policy, which balances exploration and exploitation. The MCTS algorithm explores polymorphs of boron through four steps. It begins by creating a root node from a random boron configuration. In the expansion stage, offspring nodes are generated by perturbing the structure of the root node through predefined perturbation moves,²⁸ such as changing the composition (adding/deleting atoms) and mutating the structure (swapping atom positions, deforming the structure, and straining the lattice). By conducting a finite number of random “playouts,” a qualitative measure of each node and its potential as a parent for further expansion is obtained. The playout results are then backpropagated throughout the search tree, and node selection for further expansion follows the UCB policy. The MCTS-based search is performed at 0 GPa external stress, and the objective function (enthalpy) is evaluated using the corresponding engine, LAMMPS³³ (Large-Scale Atomic/Molecular Massively Parallel Simulator), with the extended Tersoff potential³⁴ for enthalpy calculations. The pool of boron structures includes those obtained from the Materials Project³⁵ database, resulting in over 1000 structures being sampled. Degenerate structures are eliminated by comparing angular distribution functions, radial distribution functions, and SOAP³⁶ (Smooth Overlap of Atomic Positions) analysis. After a pool of structures are shortlisted based on enthalpy, they are rerun using Density Functional Theory (DFT) calculations to get more accurate enthalpy values. To compare enthalpies at high temperatures and pressures, the Gibbs free energies are estimated using a quasi-harmonic approach implemented using the phonopy³⁷ package. The methodology of the free energy calculations is implemented according to the workflow reported by Srinivasan et al.¹⁰ A detailed explanation of this method is available in the supplementary material. A pseudocode of the workflow is provided at the end of the article. The pressure–temperature grid considered for this work has six discrete pressure points (0, 5, 10, 15, 20, and 25 GPa) and a temperature grid of 0–2000 K at intervals of 10 K. Under fixed thermodynamic conditions, we determine the metastable phases by sampling the configuration space using MCTS. The ground state and metastable states, which correspond to the global and local minimum of the free energy in the configuration space, are identified by considering the lattice parameters (a, b, c, α, β, and γ) and the positions of the basis atoms (r_i). Subsequently, we compute the Gibbs free energy as a function of intensive variables (T, P) at these identified minima. The free energy and energetic ordering of the minima vary with changes in (T, P). Since the free energy of a phase (ΔG) and the extent of free energy metastability from the ground state structure (ΔG) have the same notation, we will refer to the free energy of a phase as “G” and the extent of metastability of a structure from the ground state energy as “ΔG” in this text. From the free energies (G) of each phase over the temperature range and pressure range, the free energies (G) are interpolated using a fourth order polynomial to obtain free energies at pressures steps of 1 GPa. The interpolated free energies of all the phases are compared to compute the stable phase, which is the phase with the lowest free energy (G) at a given set of temperature and pressure. Finally, we generate the complete metastable phase diagram [P′ (T, P, ΔG)], where P′ represents the most energetically favorable phase within a free energy range (ΔG) relative to the ground state at a given (T, P, ΔG). The decision boundaries (the boundaries separating the phases) are constructed using a Support Vector Machine (SVM) classifier using the scikit-learn³⁸ library in Python. More information about the decision boundary methodology is available in the supplementary material.

Figure 1 shows the schematic of the workflow¹⁰ used for the computation of the phase diagram. A pool of structures is computed via MCTS and subsequently ranked based on their enthalpies to determine the computed ground state structure and its enthalpy (ΔH_Gs). Structures are then filtered based on an enthalpy cutoff (ΔH_cutoff), which in this case is 0.2 eV, and structures having an enthalpy lesser than ΔH_Gs + ΔH_cutoff are considered as candidate structures. This cutoff was chosen because an activation barrier of roughly 0.2 eV would amount to a temperature of 2000 K (kT at 2000 K amounts to 0.17 eV). Enthalpies of the candidate structures at different pressure points are calculated along with force constants and phonon frequencies. The Gibbs free energy, of each phase at discrete (P, T) points, is approximated from the pressure–enthalpy and phonon frequency data using the quasi-harmonic approach, as shown in Fig. 2. At each (P, T) pair, the phase with the lowest Gibbs free energy is considered as the stable phase. Subsequently, for all the other phases, the metastability of that phase is calculated as the Gibbs free energy difference from the ground state structure. Decision boundaries are drawn using a SVM classifier algorithm. A 2D slice at the desired ΔG from the ground state structure gives the metastable phase diagram.

Figure 2 shows the structures of the boron phases that are considered, which are within 0.2 eV of the most stable structure. Within this cutoff, we report eight structures of boron, including two new phases, with the ground state structure being the α-phase of boron (B₁₂). The other structures in consideration are β-B₁₀₅ phase, γ-B₂₈ phase, T-B₅₀ phase, T-B₄₈ phase, R3m-B₄₅ phase, and two newly reported structures, Cm-B₅ phase and Imm2-B₇ phase. The structures reported here, barring the newly reported structures, are in good agreement with previously reported structures by computational and experimental studies.^4,39–42 α-Boron (B₁₂) exists in a rhombohedral unit cell (a = b = 4.90 Å; c = 12.56 Å). This structure consists of B₁₂ icosahedra in the corners of the rhombohedral unit cell, as shown in Fig. 2. The ground state structure reported here (α-B₁₂) is in excellent agreement with previously reported studies.^4,39–42 The next phase in consideration is β-B₁₀₅, which also crystallizes in a rhombohedral unit cell (a = b = c = 10.10 Å). It consists of clusters of boron in B₂₇ and B₁₂ icosahedra. The γ-B₂₈ phase is an orthorhombic structure (a = 5.02 Å, b = 5.59 Å, and c = 6.93 Å), consisting of B₁₂ icosahedra and B₂ dumbbells connected to each other. The T-B₅₀ phase is a tetragonal structure (a = b = 8.85 Å and c = 4.93 Å), consisting of B₁₂ icosahedra joined by a few lone boron atoms. The T-B₄₈ phase is a similar tetragonal structure (a = b = 8.64 Å and c = 4.95 Å), consisting of B₁₂ icosahedra connected to each other. Next in consideration is the R3m-B₄₅ phase, which is a rhombohedral structure (a = b = 5.81 Å; c = 11.80 Å) consisting of B₁₂ icosahedra connected to each other via B₂ dimers and B₃ triads. We report a new structure (Cm-B₅), which is a monoclinic structure (a = 5.82 Å, b = 2.90 Å, c = 4.99 Å, α = γ = 90°, and β = 55°). We report a second new structure (Imm2-B₇), which is an orthorhombic structure (a = 2.86 Å, b = 6.67 Å, and c = 5.017 Å), consisting of seven atoms. These are summarized along with their corresponding phonon dispersion curves in the supplementary material.

Figure 3 shows the relative stability curves (G vs T) of the different boron phases in consideration as a function of temperature at different pressures, namely 0, 10, 20, and 25 GPa. The curves show a convex decrease in G with increasing temperature. Under the no pressure condition, the α-, β-, and γ-phases remain extremely competitive with a difference in enthalpy in the order of 10–20 meV/atom. In the equilibrium condition at 0 K, the α-phase is the most stable phase and is subsequently replaced by the β-phase as the most stable phase at high temperatures. The extremely small difference in the enthalpies between the α-, β-, and γ-phases would explain the variation in the stability regions of the phases in the previously and currently reported regions, explained in detail later. It is also interesting to see the T-B₅₀ phase starting to become thermodynamically competitive at higher temperatures (∼2000 K). The other boron phases (Imm2-B₇, Cm-B₅, R3m-B₄₅, and T-B₄₈) remain sufficiently distanced from the stable phases at 0 GPa. One point to be noted is that the anharmonicity of the different phases is not captured in this quasi-harmonic approximation of Gibbs free energies and would lead to small variations in the stability curves and thus the phase diagrams, especially at high temperatures.

In the stability curves at 10 GPa, the α-phase and the γ-phase remain the most stable and highly competitive. However, the β-phase shows a relative increase in G with respect to the other phases, i.e., becoming more unfavorable, while the Imm2-B₇ phase shows a relative decrease in G. The Imm2-B₇ phase, R3m-B₄₅ phase, and T-B₅₀ phase become competitive among themselves. The T-B₄₈ and Cm-B₅ phases remain unfavorable. At higher pressures (P = 20 GPa), the α-phase and the γ-phase continue to remain the most stable and highly competitive. Meanwhile, the Imm2-B₇ phase becomes a little more favorable and competes with the β-phase. At high temperatures (∼2000 K), the T-B₅₀, Cm-B₅, and Imm2-B₇ phases become quite competitive, with the α-phase and γ-phase indicating an increased possibility of the formation of these phases. Finally, at pressure P = 25 GPa, we see relatively no variation in the most stable phases (α and γ). However, the Imm2-B₇ phase, at this pressure, is more favorable than the β-phase, which competes with the T-B₅₀ and Cm-B₅ phases, while the T-B₄₈ and R3m-B₄₅ phases remain unfavorable. In contrast to the stability curves at the four different pressures, we see that the α- and γ-phases remain the most stable at all pressure and temperature regimes considered. Interestingly, at each pressure, the T-B₅₀ and newly reported phases of Cm-B₅ and Imm2-B₇ phases become competitive and indicate a possibility of being synthesized. The β-phase continues to increase its unfavorability with increasing pressure, which can be chalked to the relative asymmetry of the β-phase but tends to become thermodynamically competitive at high temperatures.

From the stability curves of each phase at different pressures, the equilibrium phase diagram is obtained by compiling the P–T data of the individual phases. The decision boundaries are drawn by fitting an SVM classification model with the target classes being the unique phases found in the discrete P–T plot. Some regions of conflict are resolved manually using domain knowledge, to allow for phases to be present in a contiguous manner. The computed equilibrium phase diagram, is shown in Fig. 4(a), consisting of three phases (α-, β-, and γ-). The equilibrium phases that we report are in good agreement with the ones reported in both experimental^40,41 and computational studies,³⁹ although the pressure and temperature ranges considered in each study vary. We have chosen to limit ourselves to a temperature range of 0–2000 K, beyond which the anharmonic contributions to free energy start contributing significantly. These contributions could play the deciding role in the stability of phases, which are highly competitive, as shown in Fig. 3. Under ground state conditions (0 K and 0 GPa), our study reports that the α-phase is the most stable. This contrasts with the reported computational study,³⁹ which suggests that the β-phase is the most stable. We attribute this difference to the extremely small difference (<20 meV/atom) in the ground state enthalpies between the three phases (α, β, and γ) at 0 K, as shown in Fig. 3. The experimental study done by Parakhonskiy et al.⁴¹ depicts the α-phase as being the most stable at 1000 K and 0 GPa. Extrapolation of the equilibrium line between the α- and β-phases suggests that the α-phase would be the equilibrium ground state structure at low temperatures and pressures (∼0 K and 0 Gpa), which is in good agreement with our reported phase diagram.

At low temperatures (<1000 K), our phase diagram depicts two phases (α- and γ-) for the entire pressure range (0–25 GPa). This is in relatively good agreement with the computational study by Oganov et al.,³⁹ which shows a narrow region of the γ-phase transitioning to the α-phase and subsequently to the γ-phase. At mid-range of temperatures (1000–1500 K), we report a similar behavior as under the low temperature conditions, showing a transition from the α-phase to the γ-phase. The experimental study by Parakhonskiy et al.⁴¹ shows a transition from the β-phase to the α-phase and then to the γ-phase. We attribute this variation again to the relative low difference in the enthalpy between the α-, β-, and γ-phases. At high temperatures (1500–2000 K) and low pressures (<5 GPa), our study shows the presence of the β-phase as the most stable phase. This is in agreement with the reported computational study by Oganov et al.³⁹ and experimental studies by Parakhonskiy et al.^40,41 In this temperature range, our study shows a transition from the β-phase to the α-phase and to the γ-phase, in agreement with a reported experimental study. However, the reported studies^39–41 show a direct transition from the β-phase to the γ-phase. We attribute this discrepancy again to the relatively close enthalpies of formation of the α-, β-, and γ-phases.

Figures 4(b)–4(d) show the metastable phase diagram of boron at ΔG values of 50, 100, and 200 meV/atom, respectively, from the equilibrium ground state. Here, we reiterate that we have chosen to limit ourselves to a metastability extent of 200 meV/atom as a temperature of 2000 K would approximately result in a kinetic barrier of ∼200 meV/atom. Figure 4(b), which shows the metastable phase diagram at a ΔG of 50 meV/atom, reports the presence of the α-, β-, and γ-phases and the tetragonal T-B₅₀ phase. The tetragonal phase is present in a very narrow range at high temperatures (1700–2000 K) and low pressures (0–5 GPa). The other phases, α, β, and γ, show a drastic variation in comparison with the equilibrium phase diagram. In the low-pressure range (0–2 GPa), the closest metastable phase with respect to the ground state equilibrium phase is the γ-phase, indicating how competitive these three phases are. At relatively higher pressures, the closest metastable phase is the β-phase, which is replaced again by the γ-phase and, finally, the α-phase at high pressures (20–25 GPa). The competitiveness between these phases explains the discrepancy seen between our computed phase diagram and previously reported studies.^39–41 

Figure 4(c) shows the metastable phase diagram at a free energy difference of 100 meV/atom. The phases that are reported in this slice are the R3m-B₄₅ phase, T-B₅₀ phase, β-phase, Imm2-B₇ phase, and Cm-B₅ phase, of which only the β-phase was present in the equilibrium phase diagram. It is interesting to note that within a ΔG of just 100 meV/atom, there is a super-hard phase, T-B₅₀,⁴ suggesting that these phases are realistically synthesizable considering kinetic barriers for transformation. The newly reported phase, Imm2-B₇ phase, dominates the phase diagram in the high-pressure regime, while those of T-B₅₀, R3m-B₄₅, and β-phases occupy regions in the low-pressure regime.

At a ΔG of 200 meV/atom, the corresponding slice of the metastable phase diagram is shown in Fig. 4(d). This phase diagram shows the T-B₄₈ phase, T-B₅₀ phase, and Cm-B₅ and R3m-B₄₅ phases in equilibrium. The phase diagram is dominated mainly by the T-B₄₈ phase, leaving other regions to be occupied by the other phases, including a narrow strip at high pressures by the T-B₅₀ phase, a strip at low pressures by the Cm-B₅ phase, and a high pressure R3m-B₄₅ phase. An important point of consideration is that we have limited ourselves to phases that are 0.2 eV higher than the most stable boron phase at 0 K and 0 GPa. Consideration of more phases could potentially create regions of stability for new phases, especially at ΔG slices >150 meV/atom, because of the convex nature of free energy surfaces of the phases. However, such extensive consideration of phases is beyond the scope of the current work.

In summary, structure search using the MCTS algorithm has been performed to obtain stable and metastable boron structures. Eight structures of boron were shortlisted within an enthalpy range (ΔH_cutoff) of 0.2 eV/atom, including two metastable phases that have not been reported before. The equilibrium phase diagram is first computed, agreeing well with the existing reports. Furthermore, metastable phase diagrams of boron at ΔG of 50, 100, and 200 meV/atom are computed and regions of thermodynamic stability or metastability of these phases are explored. Our representation of metastable phases provides valuable information for experimental design and expedites the discovery of metastable phases, known for exhibiting intriguing properties that are not found in stable phases. Our framework lays the foundation for exploration and design of synthesizable metastable materials of boron.

The supplementary material accompanying the work includes the detailed descriptions on the procedure of free energy calculation, MCTS algorithm, and SVM classifier.

We thank Dr. Pierre Darancet for his contributions to this work. This work performed at the Center for Nanoscale Materials, a U.S. Department of Energy Office of Science User Facility, was supported by the U.S. DOE, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357. This material is based on the work supported by the DOE, Office of Science, BES Data, Artificial Intelligence, and Machine Learning at DOE Scientific User Facilities program (MLExchange and Digital Twins). S.K.R.S. also acknowledges the support from the UIC faculty start-up fund. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231.

The authors have no conflicts to disclose.

Karthik Balasubramanian: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Writing – original draft (equal); Software (equal). Suvo Banik: Validation (equal); Writing – review & editing (equal). Sukriti Manna: Validation (equal); Writing – review & editing (equal). Srilok Srinivasan: Software (equal). Subramanian K. R. S. Sankaranarayanan: Conceptualization (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).

The pseudocode needed to generate the metastable phase diagrams can be found at https://github.com/Karthik-Balas/elastemp-pd.⁴³