Multi-Cell Monte Carlo Method for Phase Prediction

We propose a Multi-Cell Monte Carlo algorithm, or (MC)^2, for predicting stable phases in chemically complex crystalline systems. Free atomic transfer among cells is achieved via the application of the lever rule, where an assigned molar ratio virtually controls the percentage of each cell in the overall simulation, making (MC)^2 the first successful algorithm for simulating phase coexistence in crystalline solids. During the application of this method, all energies are computed via direct Density Functional Theory calculations. We test the method by successful prediction of the stable phases of known binary systems. We then apply the method to a quaternary high entropy alloy. The method is particularly robust in predicting stable phases of multi-component systems for which phase diagrams do not exist.

phases in a chemically complex crystalline system and the composition and structure of the different phases in a single run. To the best of our knowledge, this is the first and only method that can capture the phase boundary from only one initial composition, without the need to interpolate intermediate compositions.
Coexisting supercells in Monte Carlo simulations have previously been used within the Gibbs ensemble Method of Panagiotopoulos [2] for simulation of vapor/liquid equilibrium, where atoms are constantly deleted/inserted in the cells. In crystalline solids, deleting/inserting atoms in each cell creates excessive point defects. Therefore, the Gibbs ensemble MC has not been applied to phase predictions in solids. Recently, we have introduced the first multi-cell MC relaxation method to solids, which required fixed number of atoms in each cell [3]. However, the fixed cell sizes there restrict the compositional variations and do not allow for prediction of phase fractions.
Here, we maintain the idea of a multi-cell MC but eliminate the fixed-size restriction via a modified algorithm. In the new method, each cell is assigned a molar ratio to virtually control its percentage in the total system. The molar ratios are determined by the compositions of all cells via application of the "lever rule" so that the total composition of the system is constant. In a Gibbs ensemble MC simulation, a random atom is transferred from one cell to another under constant-(µ, V, T ) condition. Here, the transfer of random atoms is achieved by changing the molar ratios of cells. Specifically, we randomly change the species of an atom in one or more cells, which we call a flip move. A flip move changes the compositions in each cell, which leads to a different set of molar ratios, or different percentages of each cell in the system. As a consequence, this is equivalent to the transfer of a group of random atoms among the cells. The total energy of each cell is computed with Density Functional Theory (DFT), with relaxation of all degrees of freedom at zero pressure. Note that the flip moves are not arbitrary and need to conserve the total number of each species among all cells.
The lever rule is often introduced in the study of binary phase diagrams [4]. It is a tool to determine the molar or volume ratio of each phase of a binary system at equilibrium at a temperature using only the compositions where the tie line crosses phase boundaries.
Assuming two supercells that represent two phases of a binary system, conservation of the initial stoichiometry requires where the superscript a indicates element a while the subscript 1 and 2 indicate the supercells/phases, n a 1 and n a 2 are the number of atoms of element a in supercell 1 and 2, respectively, n is the total number of atoms in the system and c a is the atomic concentration of element a in the system. The molar ratios (x 1 and x 2 ) are then obtained from Note that B is a constant vector, given by the initial composition of the alloy and A is updated at each MC step. Equation 1 can be easily generalized to a m-component system with m > 2.
The compositional fluctuation is achieved by changing the chemical identity of atoms rather than particle insertion/deletion. Specifically, we change the element of a randomly chosen atom to any of the elements in the system with a probability equal to the system concentration of the new element. We call these moves flips. These moves were first introduced by Kofke and Glandt to establish the semigrand canonical ensemble [5]. In our simulation, such flips can happen locally (i.e. flip a random atom in one random supercell) or globally (i.e. flip a random atom in all supercells at the same time) without causing bias. This particle flip move, combined with the application of the lever rule to multiple cells, mimics the effect of varying cell size and composition in each cell, without the need for particle insertion/deletion in the cells. This corresponds to the semigrand canonical ensemble.
On the other hand, considering all cells, the total number of atoms n as well as total numbers of each species in the entire simulation is fixed as manifested in matrix B of Equation 1. Therefore, the entire simulation represents a canonical ensemble with acceptance where ∆E tot is the change in energy of all cells combined. The contribution of each cell to the total energy is governed by the molar ratio of the cell, i.e E tot = i x i E i . Therefore, where the summation is over all the cells and primed and unprimed variables denote quantities after and before an MC step, respectively.
Accurate calculation of total energies in multicomponent alloys requires first principles calculations. We use the plane-wave-basis density functional theory code VASP [6,7], with projector-augmented-wave method pseudopotentials [8,9], the Perdew, Burke, and Ernzerhof generalized gradient approximation exchange-correlation functional [10], and a Monkhorst-Pack mesh [11] for Brillouin zone integration with Methfessel-Paxton smearing [12]. The cutoff energies used are at least 30% higher than the default values. energies obtained by the two settings at these points quantifies the error caused by using less accurate settings during the MC moves and verifies that general trends and observations remain unchanged. All the initial input structures are special quasi-random structures (SQS) [13] generated by the Alloy Theoretic Automated Toolkit (ATAT) [14] where the structures perfectly mimic those of random alloys in terms of at least first nearest neighbors.
An equiatomic Au-Pt alloy is chosen as our first benchmark, since the Au-Pt phase diagram has a miscibility gap at low temperatures [15] and AuPt alloys are thus expected to decompose into Pt and Au-enriched phases.    and 1000 K against the phase diagram by measured by Grolier et al. [16]. The points at T = 600 K (323 • C) are the average compositions of the last 10% accepted steps of each cell shown in Figure 1(b). Similarly, other points are obtained by separate (MC) 2 runs at T = 800 K and T = 1000 K. The predicted miscibility gap in this temperature range agrees very well with experimental measurements [16].
Next, we test the method with the Hf-Zr binary alloy. (MC) 2 successfully predicts the complete solubility of Hf-Zr alloy at T = 400 K, T = 700 K and T = 1000 K. Figure 3 shows study the quaternary, equiatomic, HfNbTaZr high-entropy alloy. This alloy forms a bodycentered cubic (bcc) single-phase random solid solution after casting [17].   other words these two cells do not represent a phase after the phase separation and must be discarded. Therefore, the only remaining phases are the two Nb/Ta-rich, and Hf/Zrrich ones obtained from cells 1 and 4 respectively. Since all degrees of freedom, including cell volume and shapes are optimized by VASP, the Hf/Zr-rich Cell 4 transforms into a hexagonal close packed (hcp) lattice from the initial bcc lattice. Therefore, (MC) 2 predicts the phase separation of the bcc HfZrTaNb HEA into two phases of Nb/Ta (bcc) and Hf/Zr (hcp). These predictions are consistent with those of a recent Atom Probe Tomography study on HfNbTaZr, discovering the formation of a secondary phase enriched in Hf and Zr after sufficient annealing [18].
In summary, we have introduced the (MC) 2 method to predict stable phases and phase fractions in multicomponent alloys. As we have demonstrated, each converged (MC) 2 run identifies either a region of miscibility, or the relevant phase boundaries for the simulated temperature and composition. The algorithm in (MC) 2 takes advantage of parallel computations of multiple cells and provides the unique capability of identification of relevant phases and phase boundaries without any prior knowledge of possible phases. In its current implementation, the MC algorithm takes into account the configurational entropy contributions to the free energy, while other contributions such as vibrational entropy are neglected.
This makes phase diagram predictions reliable for temperatures sufficiently below the solidus curve.