The free energy of mechanically unstable phases

Phase diagrams provide ‘roadmaps' to the possible states of matter. Their determination traditionally rests on the assumption that all phases, even unstable ones, have well-defined free energies under all conditions. However, this assumption is commonly violated in condensed phases due to mechanical instabilities. This long-standing problem impedes thermodynamic database development, as pragmatic attempts at solving this problem involve delicate extrapolations that are highly nonunique and that lack an underlying theoretical justification. Here we propose an efficient computational solution to this problem that has a simple interpretation, both as a topological partitioning of atomic configuration space and as a minimally constrained physical system. Our natural scheme smoothly extends the free energy of stable phases, without relying on extrapolation, thus enabling a formal assessment of widely used extrapolation schemes.

P hase diagrams play a central role in the physical sciences and computational methods have become a preferred route to obtain such thermodynamic information [1][2][3][4][5] . Widely used frameworks for computing phase diagrams (for example, computational thermodynamics 1,6-9 and cluster expansion [10][11][12][13][14] methods) rely on the assumption that all phases remain metastable with well-defined free energies under all conditions. Unfortunately, mechanical instabilities are common in solid-state systems [15][16][17][18] , making free energies undefined (other than by delicate extrapolation schemes). This problem has prompted a long search for reasonable and practical definitions of free energies for mechanically unstable phases 6,[16][17][18][19][20][21] and still hinders thermodynamic database development 21 . Here we propose an efficient computational solution with a natural interpretation, both as a topological partitioning of atomic configuration space and as a minimally constrained physical system. Our scheme smoothly extends the free energy of mechanically stable phases, without relying on extrapolation schemes. The results agree well with available computational and experimental estimates.
A common approach to handle mechanical instabilities is to merely assume that the functional form used to represent the data in a mechanically stable composition range will automatically provide a reasonable extrapolation into the unstable range. This 'lattice stability' picture underlies the widely popular computational thermodynamics (also known as CALculation of PHase Diagram method, or CALPHAD) and cluster expansion (also called generalized Ising model) frameworks, but is not without well-documented problems [16][17][18][19]21 . Being intrinsically noisy, extrapolations from different directions in composition-temperature-pressure space may not agree. Without careful analysis, the extrapolated free energy of an unstable phase may inadvertently lie below the free energy of a truly stable phase. Conceptually, even defining the free energy of an unstable phase is difficult. Computational approaches are being increasingly used to provide independent corroboration to the CALPHAD assessments 22 and to help address the problem of unstable phases [16][17][18][19] .
This report pursues this line of attack by proposing a general and efficient computational approach having a firmer theoretical basis that can be justified from three complementary points of view. First, it corresponds to a topological partitioning of phase space based on a simple curvature criterion. Second, it has the interpretation of stabilizing the system by constraining the minimum number degrees of freedom. Finally, the proposed definition yields free energies that vary smoothly as the system crosses the point of mechanical instability, a property that is verified both theoretically and through electronic structure calculations on actual alloy systems. The calculated free energies of unstable phases agree remarkably well with the ones obtained via existing extrapolation schemes, thus placing the latter on a stronger theoretical footing.

Results
Phase space partitioning. From a topological point of view, phase space can be naturally partitioned into configurations s based on the curvature of the system's energy hypersurface (Fig. 1). In a system of N atoms, let x denote the 3N vector of all atomic positions, let V(x) denote the potential energy of the system in that state and let c(x) be the minimum curvature at x, that is, the minimum eigenvalue of of the Hessian (the matrix of second derivatives). Hence, c(x) 40 and c(x) r0 correspond to mechanically stable and unstable regions, respectively. Within the lattice stability picture, a given lattice L is associated with many possible ways of assigning the atoms to the lattice sites. For each assignment s, we denote the ideal (or 'unrelaxed') positions of the atoms by x u s . (For each s, x u s is a 3N vector.) In the neighbourhood of each x u s , a small fraction of phase space around x u s , denoted Z s , is also associated with the configuration s, to account for static and/or dynamic displacements of the atoms around their ideal lattice positions. The free energy obtained from the classical configurational partition function associated with lattice L is then given by a 3N-dimensional integral in a familiar 'coarse-grained' form 11 : where b ¼ 1/(k B T), T is temperature and k B is Boltzmann's constant. In our scheme, the neighbourhood Z s is the largest connected set containing x u s over which the minimum curvature c(x) does not change sign. Within the set Z s , we also define x r s as the location of the minimum of V(x) within Z s (that is, the 'relaxed' positions) and Vðx r s Þ is the energy of the system in configuration s at 0 K. (Supplementary Note 1 provides more details regarding these definitions.) When a configuration s 0 corresponds to a standard mechanically stable phase, this definition agrees with the usual notion of a local potential well: if one were to 'pull' on the atoms with an increasing force, the system would break free of a given local minimum precisely at the point where a local instability develops (that is, when the minimum eigenvalue of the Hessian vanishes). For a mechanically stable phase, relaxation along a path of steepest descent from x u s 0 towards x r s 0 would always face a positive curvature along the path. In this case, the integral in equation (1) can easily be evaluated (for example, within a harmonic approximation around x r s 0 via a standard lattice dynamics calculation). When a configuration s is associated with a mechanically unstable phase, the point of minimum energy x r s in Z s is necessarily at its boundary (for otherwise the Hessian would have had to be positive definite at an interior local minimum). Finding this minimum generally involves solving a nonlinear constrained optimization problem. Checking if a given point belongs to Z s at each step of the optimization can be accomplished with the dimer method 23 , which efficiently determines the smallest eigenvalue of the Hessian. Figure 1 | Partitioning of phase space into neighbourhoods. The potential energy hypersurface V(x) (as a function of the state x of the system) defines a natural partitioning of phase space into neighbourhoods Z s , based on the sign of c(x), the local minimum curvature of V(x) (blue: negative and red: positive). Each neighbourhood (stable or not) can be assigned a welldefined free energy by integration of the partition function over that neighbourhood. In this example, the point x u s corresponds to fcc W while the point x r s 0 corresponds to bcc W (there are three symmetrically equivalent basins corresponding to bcc W) and the path joining them is the well-known Bain path. Also shown are paths of steepest descent from given 'unrelaxed' positions (x u s , x u s 0 ) towards corresponding minima (x r s , x r s 0 ) within the corresponding neighbourhoods (Z s , Z s 0 ). For the mechanically stable phase, x r s 0 lies at a local minimum while for the mechanically unstable phase, x r s lies at an inflection point along a path of steepest descent.
Fortunately, in the rather common case where the instability appears tangentially to a minimum energy path joining x u s and a neighbouring local minimum x r s 0 (Fig. 1), the point x r s of minimum energy in Z s can be easily determined as the inflection point of a path determined by the nudged elastic band method 24 . The cases where the onset of instability does not occur parallel to the path can be readily detected by performing a standard phonon analysis around the inflection point. If no other unstable mode is detected in this fashion, one can be confident that the inflection point is the appropriate point x r s (since we have thus already identified the point of instability onset). This is the situation we encountered in the cases considered in this study. Note that there may be multiple possible distinct paths, in which case the one leading to the minimum V(x r s ) determines x r s . Different possible paths can be explored by considering various specific end points of the nudged elastic band path (for example, different nearby minima). This 'path of steepest descent' picture bears some resemblance with Vineyard's transition state theory 25 , with the important distinction that an inflection point (rather than a saddle point) plays a key role in separating the different regions.
Minimally constrained system. This idea of identifying the onset of instability also leads to a second characterization of the approach as using a minimally constrained system. At x r s (and at any other point of the boundary of Z s ), the system is mechanically stable along all but one direction. (Due to point group symmetry, a few directions may simultaneously become unstable at the same point.) Since in the thermodynamic limit (N-N), the contribution of a single (or a finite number of) mode is negligible, a harmonic expansion around x r s is still possible and a standard lattice dynamics calculation can be performed (unlike the method proposed in ref. 18). (Other methods may be used to provide a higher accuracy or to handle the presence of multiple nearby minima, as detailed in Supplementary Note 1.) The thermodynamic limit thus justifies merely neglecting the few unstable modes and the system's free energy can be defined by imposing a linear constraint that freezes the few unstable modes at x r s . Our scheme can thus be viewed as constructing a 'minimally constrained' physical system. Further advantages of this interpretation are discussed in Supplementary Note 2.
In this scheme, one may be concerned that a single unstable phonon mode can be surrounded (on the same phonon branch) by infinitely many stable, but very soft, modes that may still lead to a divergence of the free energy. However, as observed in ref. 19, the free energy contribution of a phonon branch with a single point where frequency vanishes involves the integral of a singularity that is only logarithmic and thus yields a finite contribution.
Our partitioning scheme may appear somewhat oversimplified in that it only checks for the mere presence of unstable mode(s) instead of distinguishing regions that exhibit a different number of unstable modes. However, such distinction brings little more, as one has to realize that, in an extended solid, when a phonon branch touches the zero frequency axis only a small finite number of unstable modes appear. Beyond that, when a phonon branch crosses the zero frequency axis, one then has an infinite number of unstable modes. So the possible numbers of unstable modes are only 0, m, infinity, where m is a small number (that does not scale with system size). The value m is only reached at the boundaries of the regions Z s (where x r s is located for mechanically unstable phases), while the values 0 or infinity are reached inside mechanically stable or unstable regions, respectively.
Smooth extrapolation property. A third way to motivate the proposed approach is to study the behaviour of the system when some external variable a (for example, chemical composition or pressure) is continuously varied until the system crosses the point of mechanical stability at a ¼ a 0 . To this effect, we consider a parameter-dependent potential V(x,a). As shown in Supplementary Note 3, the (constrained) Helmholtz free energy of the system at positive temperature in neighbourhood Z s is a smooth function of a because it is given by the integral À b À 1 ln R x2Z s expð À bVðx; aÞÞdx in which both V(x,a) and Z s vary smoothly in a, except in rare singular cases.
The limit of zero temperature (where the free energy reduces to the energy) is useful (for example, in the context of cluster expansions) and provides the most stringent test of smoothness property (since the free energy integral reduces to the evaluation of the potential V(x) at a single point x r s , eliminating the automatic smoothing effect of the integral). To study this, we track the location of x r s as a function of a (which we denote x r (a)) and as well as the corresponding value of the potential V(x r (a),a). Under our scheme, for aoa 0 , x r (a) is at a local minimum, while for a4a 0 , x r (a) tracks an inflection point. The disappearance of the local minimum at a ¼ a 0 signals the point beyond which there no longer exists a well-defined structural energy (at 0 K) in the traditional sense. Our approach then provides a way to assign an energy via the inflection point. We can assume x to be one-dimensional without significant loss in generality by considering, once again, the potential along a path joining x u s and x r s 0 . We show in Supplementary Note 4 that V(x r (a),a) is a continuously differentiable function of a at a ¼ a 0 , so our scheme provides a smooth extension of the usual notion of local energy minimum. The fact that the minimum and the inflection point merge just indicates that the two notions agree in the limit case where they are both applicable, which is an arguably desirable property. This property is illustrated in Fig. 2 and follows the fact that, when a local minimum disappears, it does so by joining a nearby inflection point, which survives in the unstable regime (as further illustrated in Supplementary Fig. 1). This smooth extrapolation property is extremely convenient in the context of phase diagram calculations, because a series approximation (for example, by orthogonal polynomials) of the system's energy that covers both mechanically stable and unstable regions converges faster if the junction between the two regions is smooth 26 (although the junction is not sufficiently smooth to allow for a Taylor expansion across the junction). This holds for both the cluster expansion and the polynomial expansions used in CALPHAD modelling.
Application to alloy systems. We now demonstrate, through electronic structure calculations, that the above formal considerations actually lead to practically useful definitions of the free energy of mechanically unstable phases that (i) exhibit excellent smoothness properties and (ii) agree with available NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8559 ARTICLE estimates. As a benchmark, we chose the Ir-Re-W alloy system, because it is relevant to the design of substitute alloys for Re in high-temperature applications 27 and because it combines elements that each favours a different lattice: Ir is face-centered cubic (fcc), Re is hexagonal close-packed (hcp) and W is body-centered cubic (bcc). We compute the formation energies of ideal solid solutions as a function of composition, which simultaneously tests the suitability of our scheme for CALPHAD and cluster expansion modelling.
Our approach can be implemented by finding a path of steepest descent connecting the unrelaxed ideal lattice structure to the fully relaxed structure in one of two ways: (i) by using the nudged elastic band method, generalized to allow for cell shape variations 28 or (ii) by using 'damped' dynamics in which the atoms are repeatedly displaced in the direction of the force they experience by a fixed distance. The first option is useful if the relaxed structure has a lower symmetry than the unrelaxed one, while the second is computationally more efficient. In accordance with our definition of Z s and x r s , we determine the energy of the inflection point along that path, if it exists, and the energy of the fully relaxed structure, if it does not, thus suggesting the name 'inflection-detection' method. A standard phonon analysis is then performed about the point of expansion identified in the previous step. If the number of unstable modes is negligible (that is, does not scale linearly with the density of the k-point mesh), the free energy can be obtained from the phonon density of states of the stable modes in the usual way 11 . If the number of unstable modes is not negligible then a more expensive search for the approriate expansion point must be performed using the dimer method to identify points of zero minimum local curvature, from which the minimum energy point must be determined.
When searching for the nearest local minimum, one needs to allow for possible symmetry breaking (for instance, in the case of fcc W, distortions along the Bain path must be allowed so that the nearby minimum is actually bcc W). To obtain phonon entropy contributions, we simply use a harmonic approximation about x r s and a standard lattice dynamics analysis (the single unstable mode actually has null probability of falling exactly on one of the sampled k-points, so it requires no special treatment in our approach).
In Fig. 3, mechanically unstable composition ranges are readily seen where there is no overlap between the results obtained with inflection-detection and those obtained with full relaxations allowing for symmetry breaking. Formation energies obtained by relaxing the atomic geometry under symmetry constraints (such that the unrelaxed and relaxed geometries have the same space group) are seen to have poor extrapolation behaviour, because the symmetry changes discontinuously with composition (solid solutions have a low local symmetry, while pure end members have a high symmetry). Formation energies obtained under relaxations allowing for symmetry breaking exhibit a somewhat smoother behaviour, but unfortunately allow the system to move away from the intended mechanically unstable phase into the nearby mechanically stable phase (marked in parentheses). The proposed 'inflection-detection' method clearly yields the most natural extrapolation behaviour while at the same time avoiding full relaxation towards nearby mechanically stable phases. In addition, the inflection-detection predictions for pure elements agree remarkably well with the widely used Scientific Group Thermodata Europe (SGTE) values 29 , thus placing them on a stronger conceptual footing. This agreement is most likely the consequence of the smooth extrapolation property of the inflection-detection scheme. While agreement appears less good for the W hcp-bcc difference, this particular experimental estimate is based on limited data. For instance, the entropy difference between these phases is assumed to be zero in the SGTE data and an estimate which does not rely on this assumption 20 agrees significantly better.
It is instructive to further look into the case of fcc W, because it has received considerable attention. Figure 4 shows that the inflection-detection method yields values of the fcc-bcc energy and entropy difference in very good agreement with recent computations 18 . There is considerable spread in the specific CALPHAD estimates (because many postulated unstable free energies yield the same observable phase boundaries), but there is a clear 'band' of energy-entropy combinations that are generally considered consistent with observed phase boundaries and our estimate falls in the middle of this band. Supplementary Note 5 provides additional details regarding the construction of this figure and Supplementary Fig. 2 displays the same data in an alternative representation (free energy difference as a function of temperature). Over mechanically stable regions, all three calculation methods considered agree. However, the proposed inflection-detection method provides the smoothest extrapolation behaviour into mechanically unstable regions, unlike formation energies obtained with symmetry-constrained relaxations (in which the structure's space group is not allowed to change during relaxations). At the same time, the inflection-detection method avoids full relaxation towards nearby mechanically stable phases (marked in parentheses), unlike relaxations allowing for symmetry breaking (where the relaxed structures are allowed to have a lower symmetry than the corresponding unrelaxed ones). The resulting inflection-detection predictions for pure elements agree remarkably well with the widely used SGTE values 29  Discussion Caution is advised before blindly applying the proposed method to all phases exhibiting some form of mechanical instabilities. If a mechanically unstable state s can potentially relax to many multiple nearby minima, a local analysis of the free energy based on an expansion about the inflection points may need to be supplemented by an analysis of the degeneracy of the paths (which will introduce an additional entropy contribution). However, if the number of paths does not increase exponentially with system size, the corresponding entropy contribution is negligible in the thermodynamic limit. This is the case in the examples considered in this report, but would not be the case, for instance, for bcc titanium 30 . The case of bcc Ti is an example of a high-symmetry crystal structure that is mechanically unstable based on a harmonic phonon analysis but that is nevertheless 'dynamically stabilized' because the structure hops between local minima away from the high-symmetry point where mechanically instability occurs. (This can even take place in such a way that the average position of the atoms is at the high-symmetry configuration.) In such cases, the general formalism of partitioning phase space based on curvature still holds, but the dynamically stabilized phase corresponds to the collections of local mechanically stable potential wells located away from the high-symmetry point, not the mechanically unstable region containing the high-symmetry point. These many different configurations s must be accounted for to yield the correct free energy and hence a harmonic expansion about a single local minimum is insufficient. The multiple local minima can be conveniently handled via cluster expansion techniques [10][11][12][13][14] . This situation is rather different from the type of mechanical instability considered in this paper, because bcc Ti actually exists in nature while the phases we consider here (such as fcc W) do not (making the definition of their free energy more challenging). In summary, we propose a very simple and efficient method to define the free energy of mechanically unstable phases that is grounded in formal statistical mechanics and that exhibits the desirable smooth extrapolation behaviour into unstable regions that makes a direct connection with existing extrapolation schemes. Extensive testing on a particularly challenging ternary system exhibiting all three common lattices (fcc, bcc and hcp) confirms the predicted smooth behaviour and shows very good agreement with available experimental and computational estimates. The proposed method has been implemented within the Alloy Theoretic Automated Toolkit (ATAT) 31-33 as the 'robustrelax' command.
As such, the method offers promising avenue to solve the mechanical stability problem in the computational thermodynamics and cluster expansion frameworks that should benefit the large community of researchers in industry and academia that depends on phase diagrams to solve materials design problems in the automotive, aerospace, chemical process and electronic industries as well as in the general fields of materials science, solid-state physics, chemistry and planetary sciences.

Methods
Electronic structure. Electronic structure calculations were performed with the Vienna Ab initio Simulation Package (VASP) code implementing the projector augmented wave method 34 . The PBE (Perdew, Burke and Ernzerhof) functional 35 was used, with a plane wave kinetic energy cutoff of 250 eV. The k-point mesh was automatically determined, as described in ref. 31, to guarantee at least 4,000 kpoints per reciprocal atom. For ionic relaxations and force calculations (forces converged to 10 À 2 eV Å À 1 ), Fermi smearing of 0.1 eV was used with the Methfessel-Paxton scheme 36 of order 1. For total energy calculation on prerelaxed geometries, the tetrahedron method with Blöchl corrections 37 was used. Phonon calculations were performed with the supercell method using the fitfc code 33 including up to third nearest neighbour force constants in a 64 atom supercell with 0.2 Å displacements.
Solid solution model. Solid solutions were modelled using Special Quasirandom Structure (SQS) 38 , generated by the mcsqs code 39 . These SQS contain 32-48 atoms and ensure that at least the third nearest neighbour pair correlations match those of the disordered state and that all three-body correlations with a diameter below the third nearest neighbour pairs do not deviate by 41/8 from their disordered state values.
Steepest descent path determination. Nudged elastic band method calculations were performed with six images connecting an ideal unrelaxed structure to its fully relaxed counterpart. For systems with periodic boundary conditions (as considered here), the 'unrelaxed' structure is obtained by freezing all ionic degrees of freedom as well as all cell shape parameters, except for the volume, which is allowed to equilibrate (this eliminates numerical difficulties associated with very large and unphysical energy changes). To determine the fully relaxed structure, symmetry breaking was necessary for the high-symmetry pure elements. Symmetry breaking was along the Bain path for fcc-bcc instabilities and along the Burger path for bcc-hcp instabilities. For the SQS, symmetry breaking was unnecessary. When six images were found to be insufficient to resolve the energy landscape (this issue was concentrated in the bcc SQS structures), a more efficient steepest descent algorithm was used to decompose the path into 10-80 images. In all cases, piecewise cubic polynomial interpolation was used to locate the inflection point accurately.