Minimization of Gibbs Energy in High-Pressure Multiphase, Multicomponent Mixtures through Particle Swarm Optimization.

We present a global optimization method to construct phase boundaries in multicomponent mixtures by minimizing the Gibbs energy. The minimization method is, in essence, an extension of the Maxwell construction procedure that is used in single-component systems. For a given temperature, pressure, and overall mixture composition, it reveals the mole fractions of the thermodynamically stable phases and the composition of these phases. Our approach is based on particle swarm optimization (PSO), which is a gradient-free, stochastic method. It is not reliant on good initial guesses for the phase fractions and compositions, which is an important requirement for the high-pressure applications considered in this study because data on phase boundaries at high pressures tend to be extremely limited. One practical use of this method is to create equation-of-state tables needed by continuum-scale, multiphysics codes that are ubiquitous in high-pressure science. Currently, there does not exist a method to generate such tables that rigorously account for changes in phase boundaries due to mixing. We have done extensive testing to demonstrate that PSO can reliably determine the Gibbs energy minimum and can capture nontrivial features like eutectic and peritectic temperatures to produce coherent phase diagrams. As part of our testing, we have developed a PSO-based Helmholtz-energy minimization procedure that we have used to cross-check the results of the Gibbs energy minimization. We conclude with a critique of our approach and provide suggestions for future work, including a PSO-based entropy-maximization method that would enable the aforementioned continuum codes to perform on-the-fly, phase-equilibria calculations of multicomponent mixtures.

Introduction

The thermodynamics of multicomponent mixtures at high pressures is a topic of great interest in a number of fields. By components, we mean “atom types” (elements) or possibly assumed indivisible “molecular units” in certain special cases. Planetary science is awash with such systems, and prominent examples include hydrogen/helium mixtures in giant planets;[1−3] aqueous solutions of organic compounds, ammonia, carbon dioxide, and/or salts in icy moons and ocean planets;[4−8] and mixtures that may form as a result of carbon dioxide dissociation in Earth’s mantle.[9−11] Other examples are those produced from the detonation of diatomic and polyatomic species.[12−14] Moreover, there is a growing need to understand the behavior of additively manufactured alloys under high-pressure static and dynamic loading in order to certify them for use in aerospace and automobile applications.[15]

Like their single-component counterparts, multicomponent mixtures do not necessarily exist as homogeneous substances and can instead split up into two or more phases. Therefore, determining the equilibrium state for a given set of conditions requires finding which phases are thermodynamically stable and the relative amounts (mole fractions) of each phase. This procedure in single-component systems is often referred to as a Maxwell construction. The key distinguishing feature of multicomponent mixtures is that the various components may segregate preferentially into different phases, so that determining the mixture’s equilibrium state also requires finding the composition (mole fraction of every component) of each of the stable phases. This may be done by minimizing a particular free energy or maximizing the entropy, depending on the conditions that are specified. We present an optimization procedure to minimize the Gibbs energy of a multicomponent mixture, which reveals the relative amounts of the different phases and the composition of each phase, when the temperature, pressure, and overall composition of the mixture are specified. Our minimization scheme is based on particle swarm optimization (PSO), which is a stochastic method introduced in its modern-day form over 25 years ago.[16] PSO is highly versatile because it is a gradient-free method that places few requirements on the nature of the objective function that is to be minimized. As a result, it has since been applied to a wide variety of optimization problems.[17,18] A crucial feature of PSO is that it is a global optimization method which, unlike local optimization methods, does not require good initial guesses for the phase mole fractions and compositions. This is an important requirement for mixtures at high pressures because under such conditions, experimental and theoretical data on the phase diagram will likely not be available. We apply our method—which may be thought of as a multicomponent extension of a Maxwell construction—to generate phase diagrams and construct equation-of-state tables that are compatible with continuum-scale, multiphysics codes[19−21] that are widely used in high-pressure science.

We expand on the statements mentioned above by considering the case of single-component materials, for which a widely adopted thermodynamic formalism[22−25] has already been established by the high-pressure science community. In this particular formalism, the equation of state (EOS) for each phase of a single-component material is represented with a separate Helmholtz-energy model, which is subsequently divided into three contributions: a cold energy that designates the energy of the phase at 0 K, an ion-thermal free energy that accounts for the motion of the nuclei, and an electron-thermal free energy that reflects electronic excitations and ionization processes. (In some cases, it may be important to also include terms representing phenomena such as magnetism or electron–phonon interactions,[26] but this is not necessary for the vast majority of materials under most conditions.) Variants on this formalism have produced, for example, the large collection of EOSs that make up the LEOS[27−29] and SESAME[26,30,31] libraries maintained by Lawrence Livermore and Los Alamos National Laboratories, respectively. EOS construction is an iterative process that involves two distinct tasks: (1) development of phase-specific free-energy models and (2) generation of phase diagrams, the latter of which requires performing a minimization over the free-energy surface represented by the models. We focus on the second task in this study. Frequently, one is interested in a temperature–pressure diagram, but it is also common to produce a temperature–density diagram. The motivation behind making temperature–density diagrams often arises from the fact that the continuum-scale codes mentioned above read in the EOSs as two-dimensional tables of certain thermodynamic properties (such as Helmholtz energy) as a function of temperature and density, from which they can compute other properties through interpolation and numerical differentiation. These EOS tables are then used by the codes to simulate hydrodynamics, spall, and other nonequilibrium processes that may occur in materials subjected to strong disturbances like shock and ramp-compression waves.[32] Even when the kinetics of phenomena like phase transitions[33,94] or chemical reactions are sufficiently slow to preclude the establishment of local thermodynamic equilibrium, it is useful to have a capability to rigorously determine the equilibrium state because this would serve as a bounding limit for the case of infinitely fast kinetics. If the equilibrium state for a given temperature and density lies in a multiphase region, a free-energy minimization must be performed to determine the phase mole fractions, so that the total free energy of the multiphase state at that temperature and density may be tabulated. Reliable free-energy minimization methods are thus necessary for the construction of EOS tables.

In multicomponent mixtures at high pressures, the two tasks mentioned above—(1) development of free-energy models and (2) minimization over these free energies to construct phase diagrams and EOS tables—may both be sufficiently complicated that it is useful to draw a clear separation between them. High pressures naturally result in the formation of new (perhaps yet undiscovered) phases and compounds. Hence, producing a complete set of free-energy models for even a binary mixture can be extremely laborious because even before the models are to be developed, an exhaustive structure search (through experiments and/or ab initio quantum simulations) must be done to fully identify not only the fluid and terminal solid solutions, but also intermediate phases[34]—such as compounds and secondary solid solutions—that may exist over the conditions of interest. We do not intend to do that in this study for any particular mixture. Rather, our objective is to enable the second of the two tasks, namely, the construction of phase boundaries/diagrams and EOS tables via free-energy minimization. Although our focus is not necessarily on the development of free-energy models themselves, we do need such models in order to test/demonstrate our PSO-based minimization method, and for that purpose, we employ a formulation that applies well-established concepts from chemical thermodynamics[35−39] to develop mixture models that are built from the high-pressure, single-component EOSs described in the preceding paragraph. We apply this formalism to develop different example models that we subsequently use in the minimization calculations. This is somewhat akin to how one might demonstrate a new free-energy minimization technique in classical molecular dynamics simulations with representative interatomic potentials.

We organize the rest of this study in the following manner. Section and the Appendix provide a largely self-contained outline of key thermodynamic concepts and quantities relevant to mixtures that we refer to throughout the rest of the study, with an eye toward building on the very successful, high-pressure EOS framework encapsulated in the LEOS and SESAME tables. Details of our PSO algorithm and its suitability for high-pressure applications are explained in section . This is followed by section , which tests and demonstrates our PSO method with different example mixing models for the iron/gallium (Fe/Ga) system at high pressures. We show that it can reliably capture nontrivial features of the Fe/Ga phase diagrams dictated by the models, such as eutectic and peritectic points. We then demonstrate how PSO may be used to construct EOS tables of the type mentioned earlier. That is, we produce multiphase tables of Helmholtz energy as a function of temperature and density; for multicomponent mixtures, this means that each such table corresponds to a fixed overall composition. We produce tables for Fe/Ga mixtures and use them to perform EOS-related calculations (e.g., compute isentropes and Hugoniot curves), just like we may do with any table from the LEOS or SESAME libraries. Section concludes with a summary of our study. This section also presents a critical evaluation of our approach, gives suggestions on future work, and compares our approach with others that have been developed for multiphase, multicomponent mixtures.

Thermodynamic Description of Mixtures

Overview

For a multiphase mixture with different phases, its molar internal energy E, entropy S, and volume V are given bywhere χ is the mole fraction of phase i and denotes the moles of phase i divided by the moles of the mixture. In this work, we will introduce four different kinds of mole fractions, and so for the purposes of clarity, we will be very explicit when defining each of them. If phase i is at a temperature T and pressure P, we apply standard thermodynamic relations to define the molar Helmholtz energy, enthalpy, and Gibbs energy of the phase asrespectively. For the multiphase mixture as a whole, we have

If all phases are at a common temperature T and pressure P (which are fundamental requirements for the mixture to be in thermodynamic equilibrium[35]), eqs –9 lead to

To proceed further, we must have a prescription for the thermodynamic properties of the individual phases, so that we can compute G, E, S, V, F, and H of each phase i. It is useful at this point to draw a distinction between solutions and compounds. A compound is a covalently or ionically bonded substance with a fixed stoichometric ratio between its constituent elements; it remains fixed in this ratio even when it is melted to a liquid or vaporized to a gas. In contrast, a solution—which includes fluid mixtures, substitutional alloys, and interstitial alloys—may have a variable composition. If Q denotes G, E, S, V, F, or H, then Q of a compound can be represented as a linear combination of the Q of its elements (with the weights coming from the stoichiometric ratios) plus ΔQf of formation. In practice, this representation is likely to be useful only at ambient pressure, where an abundance of experimental data tends to be available even for difficult-to-compute quantities like the entropy of formation ΔSf. If we are instead interested in the compound’s properties over a wide range of pressures where such data are no longer available, a better approach is to construct a multiphase EOS specific to that compound that is separate from the EOSs for the constituent elements. Thus, if i represents a particular phase of a compound, then G, E, S, V, F, and H denote the properties extracted from the EOS corresponding to that phase. A prime example of a compound is water, for which the development of a wide-ranging, accurate EOS that covers its multiple solid and fluid phases is still an ongoing area of work. Because water has very different properties from a solution of hydrogen and oxygen, its EOS will bear little resemblance to those of its constituent elements and can be developed independently of the EOSs for those elements.

In contrast, if phase i is a solution (e.g., a fluid mixture or an alloy), a good approach is to represent G, E, S, V, F, and H in terms of the corresponding properties of the constituent components, rather than developing a separate EOS altogether like in the case of a compound. A standard prescription from chemical thermodynamics is to model the properties of each solution phase as a linear combination of the properties of the pure components that make up the phase plus additional mixing terms.[35−39] In the rest of this study, we neglect compounds and assume that each phase is a solution that is composed of up to c components. This implies that for the Fe/Ga mixture examples in section , our analysis excludes compounds like Fe3Ga, Fe6Ga5, Fe3Ga4, and FeGa3 that exist at ambient pressure[40] and are presumably still stable at higher pressures. We have noted in the Introduction that the goal of the present study is not necessarily to develop a set of free-energy models for Fe/Ga or for any other mixture. Our objective is to instead demonstrate PSO as a reliable means to perform a free-energy minimization given a set of such models. Section does this with different example Fe/Ga solution mixing models. Section includes a discussion of how the PSO algorithm may be adapted in the future to handle compounds like the ones listed above, which can then be tested once an EOS has been developed for them. Although the present study does not focus on model development per se, the formulation we outline in the rest of this section builds on existing high-pressure, single-component EOS models (e.g., those in the LEOS and SESAME tables), and so, it can be an effective basis for future development of EOSs for a diverse selection of multicomponent mixtures.

Using the index j to label the c components of phase i, we can express the molar Gibbs energy of that phase aswhere z is the mole fraction of component j in phase i and is defined as the moles of j in i divided by the total moles of phase i. Here, μ is the chemical potential of component j in phase i, and it is clear from eq why μ is also commonly referred to as the partial molar Gibbs energy of j in i. The concept of a partial molar Gibbs energy (or more generally, any partial molar property) appears prominently in the chemical thermodynamics of multicomponent solutions.[35−39] A rather general expression for μ = μ(T,P,z) isin which z = (z, z, ..., z) denotes the composition (mole fractions of the c components) of phase i, R is the gas constant, and G(T,P) is the molar Gibbs energy of component j when it exists by itself as a pure substance but is evaluated at the same T and P of the mixture. The subscript i in G indicates that this pure-component property is to be evaluated in the same state of matter as phase i. For instance, if i is a liquid, G is the molar Gibbs energy of pure j in the liquid phase.

The second and third terms on the right-hand side of eq represent contributions of mixing to the free energy. The second term, RT ln z, is the ideal-mixing contribution of j in i, while G̅excess(T,P,z) is the excess partial molar Gibbs energy of j in i and represents the nonideal-mixing contribution. It is evident that G̅excess(T,P,z) may in general be far more complicated than the ideal contribution RT ln z. For instance, the former may depend on pressure, while the latter does not. Furthermore, RT ln z depends only on the mole fraction z, while G̅excess(T,P,z) may involve the mole fractions of the other components in phase i as well. As with any EOS, developing a model for G̅excess requires fitting a set of parameters to a combination of experimental and theoretical data. For high-pressure applications, this represents a great challenge because we must fit it to data that span an enormous range of temperatures, pressures, and potentially also compositions.

Substituting eqs into 13 and 9, we obtain

The molar entropy and molar volume of each phase are given bywhere S(T,P) = −(∂G/∂T) and V(T,P) = (∂G/∂P) are the molar entropy and volume of pure j evaluated in the same state of matter (solid, liquid, etc.) as phase i, and the quantities in the square brackets are referred to as the partial molar entropy and partial molar volume of component j in phase i. The derivatives in these brackets are the excess partial molar entropy S̅excess = −(∂G̅excess/∂T) and the excess partial molar volume V̅excess = (∂G̅excess/∂P) of j in i. From these equations, one may compute the internal energy E = G + TS – PV, Helmholtz energy F = G – PV, and enthalpy H = G + TS of the individual phases. The molar internal energy E, entropy S, volume V, Helmholtz energy F, and enthalpy H of the entire multiphase, multicomponent mixture may then be calculated from eqs –3 and 7 and 8, respectively. Other thermodynamic properties of interest, such as the sound speed, may be derived by taking one or more derivatives of G, E, S, V, F, or H. This seemingly straightforward task can actually be quite complicated for the reasons explained in the Appendix, where we suggest an approximation to greatly simplify the computation of these thermodynamic derivatives.

Ideal and Nonideal Mixtures

An ideal multiphase, multicomponent mixture is defined to be one in which G̅excess = 0 for every component j and phase i over all conditions. The Gibbs energy G, Helmholtz energy F, and entropy S of such a mixture arewhile its volume V, internal energy E, and enthalpy H are

An ideal mixture is imagined as being a lattice in which the components are randomly distributed on the lattice sites. An analysis of the number of additional lattice configurations in the mixture leads to the well-known terms of the form −Rz ln z (for entropy) or RTz ln z (for free energies).[36,41] In addition, one assumes in an ideal mixture that a component experiences the same interatomic interactions, on average, in the mixture that it would experience if it were instead surrounded by other like components, that is, in the pure substance. This also implies that there is no short-range ordering (preferential mixing or bonding) or any other phenomena that lead to a deviation away from the idealized, random-distribution-on-a-lattice picture described above. As a result, V, E, and H of an ideal mixture depend only on the pure-component contributions, with no contribution from the mixing itself. [The additivity of volumes in eq has been referred to in different contexts as the “linear mixing approximation”[4,42] or as Amagat’s law.] Although the assumptions inherent in ideal mixing appear to be rather severe, if the components in the mixture are similar in nature, it can actually provide a fairly good approximation. For instance, a binary mixture of toluene and benzene (both of which are large, but similar organic molecules) behaves ideally,[35] and a recent study has found that ideal mixing may also hold for certain shock-compressed mixtures in a particular range of conditions.[42] We note that if a phase i is almost pure in a particular component j, so that z ≈ 1, then G̅excess ≈ 0 (though G̅excess for k ≠ j is not necessarily close to zero) because the environment that j experiences in this phase is very similar to what it would experience as a pure substance. In a largely phenomenological way, the goal of an accurate excess free energy model is to account for the aforementioned processes that are neglected in the ideal-mixing approximation.

We clarify that to be precise in our terminology, the definition of an ideal mixture that we have adopted is based on Raoult’s law and may be referred as an ideal Raoultian mixture.[35,38,43] In this definition, the G that appears in eq , which is the molar Gibbs energy or chemical potential of pure j, may be thought of as a reference state onto which the mixing terms are added to obtain the true chemical potential μ of j when it is mixed with the other components that make up phase i. This is the most convenient reference state to adopt for our purposes because it allows us to directly use high-pressure-oriented EOSs that have already been developed for the individual components, like the EOS tables in the LEOS and SESAME libraries, to compute G and other pure-component properties. An alternative reference state is the chemical potential of j in the ideal-gas state, for which the analogues to the excess properties are called residual properties or departure functions. (These are associated with the concepts of fugacity and fugacity coefficients, while excess properties are more closely associated with activity and activity coefficients.[35,36,44−46]) While EOSs involving departure functions are more elegant in some ways than the formulation we have chosen here, they are not convenient for our purposes because under the conditions achieved in typical high-pressure applications (e.g., ramp-compression experiments), ideal gases are often times very far away in terms of their properties from the real mixture of interest.

Application of PSO to Phase Equilibria

Algorithm Details

The main objective of this study is to develop an optimization method that enables one to construct phase diagrams for multiphase, multicomponent mixtures and produce EOS tables that can be used in high-pressure, continuum-scale simulations. In order to achieve this goal, we aim to solve the following problem: given the temperature T, pressure P, and overall composition in a mixture of phases and c components determine the mixture’s equilibrium state. The variable is the third distinct set of mole fractions that we have introduced so far; the other two sets defined earlier are the phase mole fractions and phase compositions . The mole fraction denotes the total moles of component j in the mixture (moles of j in each phase summed over all phases) divided by the total moles in the mixture. Constraining T, P, and to some fixed set of values provides a complete thermodynamic specification of the mixture, corresponding to a description of the system in the isothermal–isobaric ensemble, and therefore, it should be possible in principle to determine the mixture’s equilibrium state from this specification. The well-known condition for equilibrium in this particular situation is that the Gibbs energy G must attain a global minimum. From eq for G of an ideal mixture or the more general expression in 16 for G of a nonideal mixture, it is clear that if T and P are fixed, the only freedom one has in adjusting the value of G is to vary χ and z. Therefore, we find the equilibrium state of a mixture by adjusting the relative amounts of the different phases χ and the compositions of the phases z to minimize G subject to the constraints that all phases must remain fixed at T and P, and the following mole-balance (mass-conservation) equations

We will show later in this study that the mole-balance constraints are much easier to enforce if we introduce an auxiliary set of mole fractions , where y is defined as the moles of component j in phase i divided by the total moles of the mixture. The information content in y is the same as that in the set (χ, z), so that we can obtain the latter from the former, and vice-versa, through the following relations

The mole-balance constraints, when expressed in terms of y, take on the simpler form

Thus, the maximum value that y can take is , corresponding to the situation where all of component j resides in phase i. In each iteration of our optimization process, we convert y to χ and z, so that we can evaluate our objective function (the Gibbs energy G) according to eq , then we proceed to the next iteration by adjusting y via some prescription that obeys (29).

As stated in the Introduction, we perform the Gibbs-energy minimization using PSO, which is a stochastic method whose modern-day form is attributed to the seminal work by Kennedy and Eberhart.[16] PSO simulates the trajectory of a large number Np of particles (points) that reside in the search space. The particles together comprise the swarm. (Kennedy and Eberhart were motivated by a desire to model collective social behavior in animals, and they imagined the swarm as being a flock of birds in their original study.) In the problem of interest in this study, the search space is the Gibbs-energy surface, and each of the particles in the swarm represents—in general—a multiphase, multicomponent mixture that is fixed at the given T, P, and . If Np → ∞, the swarm signifies an isothermal–isobaric ensemble. Each particle is initially assigned to a location on the Gibbs-energy surface, meaning that it is assigned to some set of y and consequently also some set of χ and z. The y of each particle is then adjusted iteratively in a manner that depends on the particle’s current y, its best ever y (corresponding to the lowest G it has ever achieved), the best ever y achieved by any of the particles in the swarm, and potentially other quantities. The iterations are carried out until a termination criterion is met. PSO has enjoyed widespread usage since its introduction over 25 years ago, and many variants have been developed.[17,18] The assortment of PSO flavors may differ in their choice of particle initialization, the method by which the particles are moved in the search space, the termination criteria, and the various hyperparameters that influence the optimization process. A perusal of the literature reveals the perhaps obvious fact that there does not seem to be a “magic bullet”; different variants work better for different applications.

PSO has been used in a number of thermodynamic modeling studies. Cox and Christie report it to be an effective tool for optimizing EOS parameters to fit experimental data, and they have applied it to develop high-pressure-oriented, multiphase EOSs for four single-component materials: aluminum, tantalum, lead, and titanium.[47,48] PSO has also been tested along certain select points under ambient or near-ambient conditions in multicomponent systems, including mixtures of energetic materials,[49] petrochemical fluids,[50,51] and complex microemulsions.[52] Although none of these studies have applied PSO toward our goal of producing whole phase diagrams covering multiple phase boundaries of a mixture, they suggest that if extensive and careful testing is done to properly tune the optimization details, it can be a reliable method for achieving our goal.

Now that we have provided an overview of PSO, the rest of this section focuses on the specific details of our particular algorithm. The first step is to initialize the location of the particles (mixtures) along the Gibbs-energy surface by assigning mole fractions y to each. Let y denote the value of the mole fraction y for particle k at iteration l. Initializing the particle locations entails setting the mole fractions y for all phases and all components j = 1, 2, ..., c for each particle k = 1, 2, ..., Np at l = 0. To do this, we first randomize the order of the phases (so as to not bias any particular one) by creating the set , in which r is a random integer between 1 and and is different from the other integers in r. For each particle k, its mole fractions in the phase corresponding to r1 are assigned as follows:

Calculate for each component j.

Set y to be a random real number between (k – 1)dy and kdy for each j.

Thus, the k = 1 particle has mole fractions in its r1 phase that lie between 0 and dy for each component j, while the k = Np particle has mole fractions in that same phase that lie between (Np – 1)dy and for each component j. This ensures that, at least for the r1 phase, the particles in the swarm collectively provide a representative sampling of the search space. Our initialization approach is adapted from the Latin-hypercube sampling technique used in the recent study by Sterbentz et al.,[53] who have applied another stochastic method called differential evolution optimization to determine the input drive in high-pressure, ramp-compression experiments. (Incidentally, we have also tried differential evolution optimization for our problem, but we have not been able to get it to work as reliably as PSO.) If we let , so that we consider all of the remaining phases except for the phase in the sequential order given by r, we set the mole fractions of the r phase in all Np particles in a similar manner:

Calculate for each j and k.

For each j and k, set the mole fraction y to be a random real number between (q – 1)dy and qdy, where q is a random integer between 1 and Np.

Finally, the mole fractions of the phase in all Np particles are initialized by setting for each j and k. This enforces the mole-balance constraints in eq .

In every iteration, each particle’s fitness is examined by evaluating its Gibbs energy using eq . Evaluation of G requires converting each particle’s y to χ and z by applying eqs and 27. As mentioned earlier, we track y of every particle k and the best ever y achieved by any of the particles in the swarm, yswarm,best. In addition, we also track Gswarm,best, which is the lowest value of G ever achieved by any of the particles in the swarm and corresponds to yswarm,best. All of these quantities are updated, if necessary, in each iteration.

A critical part of any PSO algorithm is in updating the locations of the particles from one iteration to the next, and many of the studies in the literature[17,18,51,54] focus on proposing new ways to perform this update. If we designate the mole fractions y for particle k at iteration l as , then the mole fractions at the next iteration are given by y = y + Δy, where we calculate Δy according to the formulain which Δy from the previous iteration is set to zero for all Np particles if we are in the initialization step, where l = 0; and w, a, and b are weighting factors. All three weighting factors are based on real random numbers that lie between two specified bounds. After much testing, bounds that we have found to be effective are (−0.02, 0.05), (0, 1), and (0, 3) for w, a, and b, respectively. These random numbers are generated separately for each particle in every iteration, and so, w, a, and b will in general be different across particles and across iterations. The fact that w can be a negative number helps the algorithm escape local minima. In addition, we apply a “constriction factor”[17,18,51] on w, so that the complete expression for w iswhere lmax = 100 is the maximum number of iterations allowed (see below for a discussion on how we decided on this particular value), and rand(−0.02, 0.05) is a random number between the specified bounds. The constriction factor is always positive, and its role is to gradually decrease the magnitude of w, as PSO iterates and settles in on a minimum, with w becoming zero if l + 1 = lmax.

The mole fractions that we obtain from the sum y = y + Δy may sometimes be unphysical because one or more of the resulting y might be negative or might exceed the bounds set by the mole-balance constraints in eq . Some adjustment on y may therefore be necessary to make these mole fractions physically viable. For each particle k, we apply the following adjustment procedure on y that resembles the initialization procedure described above:

Randomize the order of the phases by creating the set as discussed earlier. Carry out the steps below for each component j.

Let , so that we traverse through all but the last phase in the order given by r and perform the following for each phase r. We remind the reader that if the phase mole fractions are set in the order dictated by r, the maximum possible value for y allowed by the mole-balance constraints is . We multiply this maximum value by scaling factors α and β—whose values we determine as explained below—to obtain  and , which we subsequently use to adjust y as necessary in the following manner:

If y < 0 (is negative), set y to a random number between 0 and .

If (exceeds the maximum), set y to a random number between and .

Finally, in order to satisfy eq , the mole fractions of the phase are determined by setting .

We have done extensive testing to determine the most effective values for the scaling factors α and β. From this testing, we conclude that it is best to set α ≈ 0 and β ≈ 1. In retrospect, this seems reasonable because if y is found to be negative, it is an indication that component j is likely not present in great amounts in phase r, and so, its value should be adjusted to be close to zero. Likewise, if y exceeds the maximum , it seems reasonable to adjust it to a number that is only slightly smaller than this maximum. Specifically, we set β = 0.99 and α = 10–9. We find that small values of α are required in single-phase regions to prevent the spurious appearance of other phases that should not be present. However, if α is too small, it reduces the efficiency in two-phase regions, especially when significant amounts of both phases are present. The value we have chosen appears to be a good compromise between these two competing needs. We have tried setting α closer to 1 and β closer to 0 and have found that this choice does not lead to lower values of Gswarm,best; in fact, this choice only reduces the efficiency of the algorithm, so that it requires more iterations to converge to the same Gswarm,best and yswarm,best.

We terminate the optimization process if one of two criteria are met. One criterion is satisfied if the maximum number of iterations lmax is exceeded. For the Fe/Ga examples that we show in section , we have found that a good compromise between reliability and computational efficiency is to set the number of particles Np = 60. With this choice of Np, we have found that PSO always terminates in less than 100 iterations, which is why we have set lmax = 100 as stated above. In single-phase regions, it usually converges within 15–20 iterations, while in two-phase regions, it usually does so in under 40 iterations, but in rare instances can take up to 80 iterations. A second termination criterion we check is if the method stops producing meaningful improvements. We define this as having occurred when Gswarm,best at the current iteration minus the running average of this quantity over the past 10 iterations converges to within a very tight tolerance of 10–10RT. We evaluate whether or not this criterion is met starting in the 15th iteration.

Once the optimization process has terminated, we use eqs and 27 to translate yswarm,best to χ and z, and we take the resulting set (χ, z) as being the equilibrium state that corresponds to the given T, P, and . If component j is insoluble in phase i, this will be revealed by the fact that the equilibrium value of the mole fraction z will be a very small number, typically 10–9 or less. Likewise, if phase i is not present in the equilibrium state (this state may lie in a single-phase region for instance), χ will be a negligibly small number. From knowledge of the equilibrium χ and z, one can compute other properties of the mixture (internal energy, density, sound speed, etc.) by applying the formulas in section and the Appendix to the chosen EOS for the mixture of interest.

Suitability of PSO for High-Pressure Phase Equilibria

A key set of equations that we use in our testing of PSO is the equality of chemical potentials in eq , which is a fundamental thermodynamic requirement for equilibrium.[35−39] Suppose that the equilibrium state for a given is an -phase state, where is an integer that is less than or equal to the maximum possible number of phases . The chemical potentials across the various phases must then satisfy

We may solve these equations with a direct numerical method (e.g., a Newton–Raphson solve) to obtain the phase compositions . A crucial point we emphasize is that we do not explicitly enforce these constraints in PSO, but the equilibrium state that it converges to must nevertheless satisfy these equations because they are necessary but not sufficient conditions for the Gibbs energy to achieve a minimum. Thus, they serve as an important check on the phase boundaries predicted by PSO.

Direct numerical solutions of eq offer a potentially simple and computationally efficient means to infer phase boundaries and construct phase diagrams, but they have certain limitations that a more sophisticated approach like PSO is able to overcome. One limitation is that the direct solutions reveal information only about the phase compositions z, but not the relative amounts of the different phases encapsulated in the mole fractions χ, whereas PSO reveals the complete set (χ, z). Another issue is that a direct method for solving eq will converge to a sensible solution only if there actually are phases in equilibrium at the given . It will fail to converge if this is not the case, especially if a gradient-based method is used. This means that number of phases must be known a priori if using a direct method. Contrast this with PSO, which requires input only about the maximum possible number of phases and free-energy models for these phases. Given this input, PSO will determine the actual number of phases in the equilibrium state and will thereby eliminate much of the guesswork that would be required if one were to instead piece together a mixture’s phase diagram from the direct solutions alone.

Furthermore, it is important to reiterate that the equations in eq represent necessary but not sufficient conditions for equilibrium; they cannot distinguish between local minima and the global minimum. Direct solutions of those equations therefore cannot resolve ambiguities when there are multiple different sets of phases that can be in equilibrium (each set represents a different minimum) at the given T, P, and , meaning multiple different ways for the resulting phase boundaries to combine together to form the phase diagram. An example of this kind of ambiguity will be presented in Figure . The direct solutions tend to be useful only when the “right answer” is known in some sense. For instance, we will apply them in this study as an a posteriori check on the PSO-predicted phase diagrams. Direct solutions also enjoy usage when experimental data on the phase diagram are readily available, which is the case for some mixtures at ambient pressure. In such situations, eq may be solved directly to verify whether phase boundaries predicted by EOS models developed for these mixtures can reproduce the experimental data. One of the authors of this paper has in fact done this[55,56] for condensed-phase mixtures of unreacted high explosives[57] at ambient and near-ambient pressures. However, in the Introduction, we have alluded to a quite different reality that exists at higher pressures, where data regarding phase boundaries are not available for the vast majority of mixtures. For these mixtures, the only such information that one may reasonably hope to obtain are predictions from atomistic simulations at a few select points along a particular phase boundary. Consequently, the free-energy models for these mixtures will have to be fit to more readily available data, such as isotherms, Hugoniot curves, and isentropes. Thus for high-pressure applications, one must be able to reliably generate phase diagrams with little a priori knowledge about the appearance of the phase diagram, that is, how the phase boundaries are connected together. This is precisely what a global optimization method like PSO, which does not require one to specify good guesses for the phase fractions χ and compositions z, is intended to do. As long as it is given EOS models for the different phases of the mixture that might be present, it will automatically sort out the different possibilities for the phase boundaries to find the one that minimizes the Gibbs energy represented by those models.

Figure 7

Gibbs energy relative to the PSO-predicted results as a function of composition for Fe/Ga mixtures at different fixed temperatures and pressures: (a) 4500 K and 205 GPa; (b) same as in (a), except that the y-axis is on a linear scale; (c) 3450 K and 205 GPa; and (d) 7150 K and 700 GPa. The orange curves indicate the Gibbs energy of the fluid phase minus that of the PSO result, normalized by RT. The blue and green curves show the same ΔG/RT comparison for the Ga-III/Fe and Ga-IV/Fe phases, respectively. The black curves portray the Gibbs-energy difference between the various hypothetical scenarios described in the text and the PSO result. These scenarios access local minima along the Gibbs-energy surface. The vertical dashed lines represent boundaries between two-phase (fluid–solid or solid–solid) regions for the PSO (equilibrium/global minimum) result, while the vertical dotted lines show the same for the hypothetical scenarios. The dashed–dotted vertical lines in (d) show the left (hypoeutectic) branch of the fluid–solid phase boundaries in Figure ; these boundaries are traversed in both the PSO and hypothetical scenarios.

Figure 2

Temperature–composition phase diagrams of ideal Fe/Ga mixtures in which the pressure P is fixed at some value in each figure: (a) P = 205 GPa; (b) magnified view of (a) in the Ga-rich side; (c) P = 700 GPa; and (d) magnified view of (c) in the Ga-rich side. The x-axis refers to the overall mole fraction of Fe, which we symbolize as , so that the overall mole fraction of Ga is . Each “+” symbol therefore denotes a particular set for which we have carried out a Gibbs-energy minimization using PSO. The Gibbs energy in all cases here is computed from eq . We consider three possible phases for the Fe/Ga mixtures: (1) fluid (orange), (2) Ga-III/Fe (blue), and (3) Ga-IV/Fe (green). States that lie in the two-phase, fluid–solid region are colored in red. The liquidus and solidus curves that enclose this lens-shaped region are illustrated by the black solid curves. These curves indicate the composition z of the fluid and solid phases. The phase mole fractions χ in each of the two-phase states are computed directly by PSO but can also be determined graphically from the liquidus and solidus curves by applying the lever rule. The horizontal dotted lines in (a) and (b) reflect the transition temperature between pure Ga-III and Ga-IV at 205 GPa.

Figure 3

Temperature–composition phase diagrams of Fe/Ga mixtures at a fixed pressure of 205 GPa: (a) the entire composition range; (b) magnified view of the Fe-rich side; and (c) magnified view of the Ga-rich side. We have used the same coloring scheme for the different phases/regions as in Figure . Like in that figure, the two solid phases (Ga-III/Fe and Ga-IV/Fe) are modeled as ideal mixtures, but unlike in that figure, the fluid phase is not an ideal mixture but instead obeys the nonideal-mixing model in eqs and 44. At the 205 GPa pressure represented in this figure, the resulting excess Gibbs energy in the fluid is positive, meaning that nonideal mixing destabilizes the fluid, making it less thermodynamically favorable. The boundaries of the fluid–solid region (the liquidus and solidus curves) are indicated by the solid curves. For comparison, we have superimposed these same boundaries from the ideal-mixing case of Figure a,b (see the dashed curves). The horizontal dotted lines in (a,c) again reflect the Ga-III/Ga-IV transition temperature at 205 GPa.

Figure 5

Temperature–composition phase diagrams of Fe/Ga mixtures at a fixed pressure of 205 GPa: (a) the entire composition range; (b) magnified view of the Fe-rich side; (c) magnified view of the Ga-rich side; and (d) even more magnified view deep in the Ga-rich side, near the peritectic point. The Ga-III/Fe phase is treated as being ideal, while the fluid phase is described with the same excess Gibbs energy model as in Figure . The difference with that figure is that we now model the Ga-IV/Fe phase with the nonideal-mixing free energy described in the text, which stabilizes that phase. (Contrast this with the situation in the fluid, where nonideal mixing destabilizes it at 205 GPa.) The Ga-III/Fe and Ga-IV/Fe phases are separated by a two-phase, solid–solid region that we have colored in cyan. The solvus curves that represent the boundaries of this region are depicted by the dashed lines. This particular phase diagram has a peritectic point at a temperature and Fe overall mole fraction of about 3458 K and 0.0083, respectively. This point represents a three-phase state, where Ga-III/Fe at that composition is in equilibrium with both the fluid and Ga-IV/Fe. The tie line connecting the three phases at the peritectic temperature is illustrated by the dotted line in (d).

Figure 6

Temperature–composition phase diagrams of Fe/Ga mixtures at a fixed pressure of 700 GPa: (a) the entire composition range; (b) magnified view of the Fe-rich side; and (c) magnified view of the Ga-rich side, near the eutectic point. At this pressure, the nonideal-mixing model for Ga-IV/Fe (the same as the one we have used in Figure ) destabilizes that phase, so that it is displaced almost completely by the Ga-III/Fe phase. Like in previous figures, the liquidus and solidus curves that enclose the fluid–solid regions are represented by the solid curves, while the solvus curves that form the boundaries of the solid–solid region are illustrated by the dashed lines. This particular phase diagram has a eutectic point at a temperature and of about 7106 K and 0.188, respectively. This point represents a three-phase state, where the fluid phase at that composition is in equilibrium with both Ga-IV/Fe (which has a slightly less Fe-rich composition) and Ga-III/Fe (which has a more Fe-rich composition).

Results and Discussion

Here, we demonstrate our PSO algorithm on the Fe/Ga system, using LEOS tables for pure Fe and pure Ga combined with different example models for the mixing free energy. We verify through different ways that it can reliably locate the minimum on the Gibbs-energy surface fashioned by these models. We illustrate this by producing temperature–composition phase diagrams at different fixed pressures. (Although not shown in this study, we could also apply the same algorithm to fix the temperature and generate pressure–composition diagrams instead.) Finally, we explain how PSO may be used to construct multiphase LEOS- or SESAME-style tables for a mixture at a particular overall composition , again using the Fe/Ga system to demonstrate the construction process. In this process, we first produce a temperature–pressure phase diagram for the mixture using PSO and then apply post-processing steps to obtain a table of Helmholtz energy versus temperature and density, which may be graphically visualized in the form of a temperature–density diagram. We compare Hugoniot and isentrope curves obtained from the resulting Fe/Ga tables with those from the pure Fe table.

Before proceeding to the mixtures, we briefly introduce the EOSs we use for pure Fe and Ga, along with their phase diagrams. For Fe, we use a new two-phase EOS developed by Benedict et al.,[58] which is intended for temperatures above roughly 3200 K and pressures above around 100 GPa. A recent study on shock-compressed Fe reports that the solid under these conditions is in the hexagonal closed-pack (hcp) phase.[59] As a result, the two phases modeled by the EOS developed by Benedict et al. are fluid and hcp. To keep the nomenclature simple, we use “fluid” as a general term to indicate liquid or supercritical fluid or any other fluid-like phase, all of which are represented by the same set of free-energy models in the EOS. The temperature–pressure phase diagram predicted by this EOS is shown in Figure a. For Ga, we use the five-phase EOS developed by Wu et al.,[60] though only three phases—the fluid plus two solids, Ga-III (which has a body-centered tetragonal [bct] structure) and Ga-IV (face-centered cubic [fcc])—are thermodynamically stable in the applicable temperature–pressure range of the Fe EOS developed by Benedict et al. The phase diagram of the Wu et al. Ga EOS is shown in Figure b–d.

Figure 1

Temperature–pressure phase diagrams of iron (Fe) and gallium (Ga) constructed from the pure-component EOSs used in this study: (a) Fe; (b) Ga; (c) same as (b), except that here we have neglected the fluid phase, so that metastable extensions of the solid phases can be seen; and (d) magnified view of (c) near the lower end of the temperature and pressure range. Each “+” symbol denotes a point at which we have performed a free-energy minimization to determine the thermodynamically stable phase(s). The Fe phase diagram is produced from the two-phase EOS developed by Benedict et al.,[58] where we have used blue and orange to indicate the hcp and fluid phases, respectively. We model Ga with the five-phase EOS developed by Wu et al.,[60] though only three phases (Ga-III = blue, Ga-IV = green, and fluid = orange) are stable under the conditions in this figure. The fluid–solid (melt) phase boundary in all figures is illustrated by the dashed curves, while the dotted curves in (c) and (d) indicate the Ga-III/Ga-IV phase boundary.

Temperature–Composition Phase Diagrams

Ideal Mixtures

When considering a mixture of two or more components, a question naturally arises as to what are the possible number of phases and the identity of these phases. In the Fe/Ga system at the pressures shown in Figure , obvious choices are the fluid phase, which is a mixture of fluid Fe and fluid Ga, and three terminal solid solutions: hcp, bct, and fcc. In the first of these solid solutions, Fe in the hcp phase would serve as a solvent for Ga solutes in that same phase. Likewise, the bct and fcc phases would be composed of Fe solute in those phases dissolved in a Ga-III or Ga-IV solvent, respectively. At intermediate compositions, there might also be compounds (e.g., Fe3Ga, Fe6Ga5, Fe3Ga4, and FeGa3) and potential secondary solid solutions—such as a body-centered cubic (bcc) phase—that have yet to be discovered. Developing a complete set of free-energy models to cover all of these possibilities is not the purpose of this study, nor is it necessary to demonstrate the efficacy of PSO. Given the lack of currently available EOSs (e.g., high-pressure EOSs have not yet been developed for any of the compounds, nor do we have a bcc-phase EOS for either Fe or Ga that is usable under such conditions), we take in which the three phases include the fluid and two solid solutions that we refer to as the Ga-III/Fe and Ga-IV/Fe phases. In both of these solid solutions, the properties of pure Fe are modeled with the hcp-phase EOS from Benedict et al. We present various mixing models and demonstrate that PSO may be used to reliably locate the minimum on the Gibbs-energy surface represented by these models to produce coherent phase diagrams.

The simplest test of PSO is the case of ideal mixing. Under this assumption, the Gibbs energy that we minimize is taken from eq and may be expressed aswhere the phase index i covers the fluid, Ga-III/Fe, and Ga-IV/Fe phases (in that order), while the component index j = 1 stands for Fe and j = 2 denotes Ga. We have used PSO to produce temperature–composition diagrams at two different pressures: 205 GPa and 700 GPa; see Figure . The fluid–solid region (colored in red) in this figure is lens-shaped, which is a characteristic feature of ideal mixing. In these two-phase regions, the equilibrium conditions in eq require that the chemical potentials satisfyin Figure a,b, where the fluid (i = 1) is in equilibrium with Ga-III/Fe (i = 2), while they must satisfyin Figure c,d, where the fluid is in equilibrium with Ga-IV/Fe (i = 3). The PSO-predicted two-phase regions in Figure satisfy these equations, without requiring any explicit enforcement, to within 10–5RT (the left- and right-hand sides agree to about 0.00005% or better). We have also verified that the liquidus and solidus curves illustrated in the figure are essentially indistinguishable from what we would obtain if eq or eq were explicitly solved with a direct numerical method to a very tight tolerance of 10–9RT.

PSO is also able to correctly describe the solid-phase equilibria dictated by the Gibbs-energy surface in eq . We can see from Figure d that at 205 GPa, Ga-IV is more thermodynamically stable than Ga-III at low temperatures, but there is a relatively narrow temperature band of about 100 K directly below the melt temperature, where Ga-III is more stable than Ga-IV. Thus, we would expect both the Ga-III/Fe and Ga-IV/Fe phases to be present in the temperature–composition diagram at this pressure, and that is indeed what we observe in Figure a,b. These PSO-generated phase diagrams suggest that there is no two-phase, solid–solid region, and instead, there is an abrupt transition from Ga-IV/Fe to Ga-III/Fe or vice-versa. The absence of a solid–solid region may be understood from the following arguments. If such a region were to exist, all points in it would satisfywhich means that for Fe (j = 1), we have from eq thatand because z22 = 1 – z21 and z32 = 1 – z31, we have for Ga (j = 2) that

Because the pure-component Gibbs energies of Fe in both solid solutions are described by the same free-energy model (see our discussion at the beginning of this subsection), we have G21 = G31, which therefore implies that z21 = z31 and G22 = G32. The last equality indicates that under the assumption of ideal mixing, the two solid solutions can be equilibrium only at the transition temperature between pure Ga-III and Ga-IV, which is represented by the dotted lines in Figure a,b. Moreover, at 700 GPa, the Ga-III/Fe phase is absent, as we would expect from Figure .

The phase diagrams that we generate in this study are produced through an iterative procedure that resembles the process of how one might put together a jigsaw puzzle. In the first iteration, a set of overall mole fractions is chosen, and for each , we perform PSO calculations over a coarsely spaced set of temperatures. This produces a rough overview of the phase diagram. In subsequent iterations, we refine it by performing calculations over additional choices of and/or homing in on multiphase regions with a more closely spaced temperature grid. Two iterations are sufficient for the phase diagrams in Figure , but multiple rounds of refinement may be needed for more complicated ones. Results for some of the points are omitted in order to improve the visual appearance/clarity.

Nonideal Mixtures

If we relax the assumption of ideal mixing, the expression for the Gibbs energy G becomes

It is conceptually useful to decompose 39 so thatin which ΔGpure = ΔGpure(T,P,z) is the change in Gibbs energy that occurs when the components of phase i are physically brought together in the proportions indicated by zand ΔGmix = ΔGmix(T,P,z) is the change in Gibbs energy that occurs when the components are chemically mixed together (dissolved) in phase i

In order for any functional form of ΔGmix to be physically valid, it must go to zero if phase i is pure in any component j, so that z = 1 and z = 0 for k ≠ j. Furthermore, the nonideal-mixing effects encapsulated in the excess Gibbs energy must vanish in the limit, as T → ∞. The excess Gibbs energy may be nonzero at T = 0, reflecting “cold” contributions to the mixing that may be computed from ab initio quantum simulations or inferred from high-pressure cell experiments conducted at low temperatures. For the Fe/Ga examples in this study, we employ the functional formHere, nonideal mixing is represented in terms of a Margules model[35,43] of the form GMzz, and we set the temperature and pressure dependence of the Margules coefficient GM aswhere η, σ, and θ are phase-specific constants.

As a first step toward nonideal mixing, we consider the situation where the fluid phase is a nonideal mixture described by eqs and 44, while the two solid solutions still undergo ideal mixing (which essentially implies that for the time being, we set η2 = η3 = 0). We set η1 = 0.5 cm3/mol, σ1 = 100 GPa, and θ1 = 2000 K for the fluid. With this choice of σ1 and η1, the Margules coefficient GM is positive (negative) for all temperatures and compositions if the pressure P is below (above) 100 exp(1) ≈ 271.8 GPa. Therefore, we expect nonideal mixing to destabilize (and thus contract the stability field of) the fluid at 205 GPa, while it should have the opposite effect at 700 GPa. This is exactly the behavior depicted in Figures and 4, respectively. Again, we can validate the results predicted by PSO with a direct numerical solution of the equilibrium conditions in eq or 35. The liquidus and solidus curves that we obtain from this direct solution overlap with those predicted by PSO. This includes the behavior on both sides of the minimum at 700 GPa that occurs around an Fe mole fraction of , which is a particular two-phase state where the coexisting fluid and solid [Ga-IV/Fe] happen to have the same composition.

Figure 4

Temperature–composition phase diagrams of Fe/Ga mixtures at a fixed pressure of 700 GPa: (a) the entire composition range; (b) magnified view of the Fe-rich side; and (c) magnified view of the Ga-rich side. Here, we have applied the same models and assumptions as in Figure , with the key difference being that now at 700 GPa, our nonideal-mixing model for the fluid stabilizes that phase and thereby expands its stability field. Again, the solid curves show the fluid–solid phase boundaries in the nonideal case, while the dashed curves are taken from Figure c,d and portray the same boundaries if we instead assume ideal mixing.

We may add a degree of complexity by allowing for nonideality in one of the two solid solutions, the Ga-IV/Fe phase. We model the Gibbs energy of mixing in Ga-IV/Fe with eqs and 44, and we set two out of the three parameter values to be the same as in the fluid, so that σ3 = σ1 = 100 GPa and θ3 = θ1 = 2000 K. However, we choose η3 = −0.1 cm3/mol, and because η3 has a different sign and is smaller in magnitude than η1 = 0.5 cm3/mol, this particular nonideal-mixing model has the opposite, though weaker, effect as in the fluid. Therefore, nonideal mixing at 205 GPa increases the thermodynamic stability of Ga-IV/Fe, while it destabilizes Ga-IV/Fe at 700 GPa. Results at these two pressures are presented in Figures and 6, respectively. Comparing Figures with 3, both of which are at 205 GPa, we see that nonideal mixing in Ga-IV/Fe at this pressure results in a significant expansion in the stability field of that phase, to the extent that it almost entirely displaces Ga-III/Fe. The Ga-III/Fe phase appears in only a narrow range of low Fe concentrations (high Ga concentrations), as shown in Figure d. Ga-III/Fe and Ga-IV/Fe are separated by a two-phase, solid–solid region that we have colored in cyan. All points in this region satisfy the equality of chemical potentials in eq , which when we substitute the Margules model in eqs and 44 for the excess Gibbs energy, may be written asand

The behavior of the solid-phase equilibria in Figure may be rationalized by considering the term GMz312 that appears in eq , which is negative at 205 GPa (because GM < 0 at this pressure) for all compositions. At very low Fe concentrations where z31 ≈ 0, this term is small in magnitude, and so, it cannot compensate for the fact that the Gibbs energy G32 of pure Ga-IV is larger than G22 of pure Ga-III for most temperatures at this pressure. This explains why Ga-III/Fe persists at very low Fe concentrations despite the stabilizing effect of nonideal mixing on Ga-IV/Fe. The opposite occurs at higher Fe concentrations, where GMz312 becomes too large in magnitude, such that Ga-IV/Fe is favored over Ga-III/Fe for all temperatures shown in the phase diagram. The solid–solid region lies in an intermediate range of compositions and represents states for which it is possible to simultaneously satisfy both eq and eq . The intersection of the solid–solid and fluid–solid regions results in a peritectic point, which is a state where all three phases are in equilibrium. At temperatures above (below) the peritectic temperature [the dotted tie line in Figure d], the solid that is present along the solidus is the Ga-IV/Fe (Ga-III/Fe) phase.

Similar reasoning may be applied toward comparing Figures and 4, both of which are at 700 GPa. Nonideal mixing at this pressure destabilizes Ga-IV/Fe, so that it gets displaced by Ga-III/Fe throughout most of the composition range. Ga-IV/Fe is present only at low Fe concentrations, where the magnitude of the excess mixing term GMz312 that appears in eq is small. The liquidus curve in Figure exhibits a minimum, and when combined with the two-phase, solid–solid region, the result is a eutectic point that can be seen in Figure c. This is a three-phase equilibrium state that represents lowest temperature where the fluid can exist at 700 GPa.

Like we have done for Figures –4, we can compare the PSO-predicted phase diagrams in Figures and 6 against direct numerical solutions of the equilibrium conditions in eqs –36. Solving those equations produces the liquidus, solidus, and solvus curves that enclose the two-phase regions, and we have verified that in all cases, results for these phase boundaries from the two methods are indistinguishable from each other. However, because Figures and 6 have a fair degree of complexity, and there are different possibilities for these phase diagrams, it is worth doing a more thorough investigation to ensure that PSO does indeed yield the true equilibrium state. For example, Figure shows that above the peritectic temperature, Ga-IV/Fe dominates over Ga-III/Fe due to the particular nonideal-mixing model we have chosen. If we assume that the phase along the solidus is actually Ga-III/Fe instead of Ga-IV/Fe, the phase boundaries that we obtain from a direct numerical solution—which represent local minima in the Gibbs energy—turn out to be similar (displaced toward higher Fe concentrations by about 2 mol % or less) to those displayed in the figure. Perhaps surprisingly, the Gibbs energies along the liquidus and solidus in that hypothetical scenario are actually lower for most conditions than the Gibbs energies along those same boundaries in the figure. Thus, based on the direct numerical solution alone, one might be led to conclude that Ga-III/Fe is actually the true equilibrium phase along the solidus. However, an examination of the Gibbs energy as a function of composition at a fixed temperature reveals that this is not the case. Figure a,b illustrates results for the 4500 K isotherm in Figure . The hypothetical scenario that we have laid out above, which is represented by the black solid curve, is always equal to or higher in Gibbs energy than the PSO result. This occurs because the liquidus and solidus in the PSO result (the vertical dashed lines) are shifted to the left compared to those boundaries in the hypothetical scenario (the vertical dotted lines). As a result, for any set of that lies within both of these two-phase regions, the solid phase fraction will be larger in the PSO case than in the hypothetical scenario. The differing relative amounts of fluid versus solid works out such that the two-phase state predicted by PSO has a lower Gibbs energy than the two-phase state in the hypothetical scenario, despite the fact that the latter may have lower free energies along its liquidus and solidus. Figure c shows a similar analysis for 3450 K and 205 GPa, which is slightly below the peritectic temperature of 3458 K at that pressure. The black curve in this figure represents the hypothetical scenario where we ignore the solid–solid transition and assume that the phase that is present along the solidus is Ga-IV/Fe instead of Ga-III/Fe. Figure d demonstrates results for the 7150 K isotherm at 700 GPa, which lies a little above the eutectic temperature of 7106 K at that pressure (see Figure ). Here, the black curve denotes the local minimum, where the branch of the solidus that lies to the right of the eutectic point is composed of Ga-IV/Fe rather than Ga-III/Fe. We note that in all the hypothetical scenarios depicted in Figure , the free-energy differences are rather small; even well into the two-phase regions, ΔG/RT for the black curves in most instances is less than 10–3RT. This suggests that PSO reliably yields the true equilibrium state even when there are competing local minima that are close in free energy to the equilibrium state. It illustrates the point we mentioned in the previous section about how PSO, being a global optimization scheme, is able to resolve ambiguities in phase boundaries and diagrams, whereas direct numerical solutions of eqs –36 [and (32) in the general case] cannot do this.

As an additional way to verify the results of the Gibbs-energy minimization, we have developed a separate version of PSO that minimizes the Helmholtz energy instead and have used it to cross-check the results of Gibbs-energy minimization. The details and motivation behind this second PSO method is best explained in the context of EOS table generation, which is the subject to which we now proceed.

Generation of EOS Tables

We have discussed in the Introduction the practical need for EOS tables in high-pressure science. An EOS table for a single-component system is simply an array of numbers representing the Helmholtz energy of the equilibrium (possibly multiphase) state on a grid of temperatures T and densities ρ. For a mixture, we would also have to fix the overall composition . Thus, table construction involves finding the equilibrium state for a given set of , much like how we have used PSO to find the equilibrium state for a given set of . Basic thermodynamics dictates that one must minimize the Helmholtz energy in the former situation. Currently, there does not exist a method to rigorously produce multiphase EOS tables for mixtures in this manner, and so, we have explored PSO as a means to do that. Our algorithm for Helmholtz-energy minimization is very similar to that for Gibbs-energy minimization outlined in section (indeed, this flexibility of PSO is one of the reasons we have chosen it in the first place), with the swarm representing a canonical ensemble instead of an isothermal–isobaric ensemble. The only procedural difference is that Helmholtz-energy minimization requires one additional step, where the pressure P of each particle in the swarm—here, each particle represents a multiphase, multicomponent mixture that is fixed at the specified T, ρ, and —must be computed in every iteration. This P is common to all phases in the particle and is the pressure corresponding to the particle’s χ and z at the current iteration plus also the specified T and ρ. Unfortunately, while this additional step is conceptually simple, it is computationally expensive. A minimization over one set of for the Fe/Ga examples is 2 orders of magnitude slower than an equivalent calculation for a set in Gibbs-energy minimization. As a result, our Helmholtz-energy minimization algorithm cannot, in its present form, construct multiphase EOS tables for mixtures in a reasonable amount of time. This may become possible with the speedup options discussed in section , but its usage is currently limited to serving as a cross-check for the Gibbs-energy minimization through the following process:

Select a particular and perform a Gibbs-energy minimization with PSO to obtain the equilibrium values of the phase fractions χ and compositions z.

Input this χ and z, as well as T and P, into the given free-energy models to compute ρ.

Perform a Helmholtz-energy minimization with PSO on and verify that its predicted equilibrium state reproduces the P, χ, and z from step 1.

We have cross-checked several points with both ideal and nonideal-mixing models and have found that the two minimization procedures yield the same equilibrium state for all the points that we have tested.

It is still possible to produce EOS tables with the more indirect procedure that we now describe. In this procedure, we first produce a temperature–pressure phase diagram for a fixed overall composition through Gibbs-energy minimization with PSO. Figure shows such diagrams for the case of a 90% Fe and 10% Ga mixture. A feature we point out is that the one-dimensional melt curve in a single-component system like pure Fe transforms to a two-dimensional, fluid–solid region in the binary Fe/Ga mixture. (These two-phase regions in iron–silicate alloys have practical importance in studying the structure of the inner Earth.[61]) The increased dimensionality compared to the single-component case is a natural consequence of mixing and is manifested in the well-known Gibbs phase rule.[35−39,95] Simple mixing protocols that combine different single-component EOSs without properly accounting for changes in phase boundaries would inherently not be able to capture the increased dimensionality. After the temperature–pressure diagrams are produced, we use our free-energy models to compute different properties like the density ρ and Helmholtz energy F at each T and P point, and we fit these results to two-dimensional splines of T and P. We find the lowest and highest densities spanned by these diagrams, which are represented by the dashed curves in Figure . We then employ the splines over the temperature and density range to generate a table of F as a function of T and ρ for the specified . This table will necessarily be multiphase and will properly reflect the increased dimensionality (in the Gibbs phase rule sense) that arises from mixing.

Figure 8

Temperature–pressure phase diagrams for a fixed overall composition of 90% Fe and 10% Ga where (a) all three solution phases form ideal mixtures and (b) the fluid and Ga-IV/Fe phases are described by nonideal-mixing models in the previous section, while Ga-III/Fe undergoes ideal mixing. The various regions are labeled with the same coloring scheme used in all of the previous phase diagrams. We follow the procedure explained in the text to convert these temperature–pressure diagrams to multiphase EOS tables of Helmholtz energy vs temperature and density. An additional round of PSO calculations could be performed to more clearly refine the boundaries of the fluid–solid region colored in red, but doing so is not necessary to demonstrate the table-construction process. The left and right dashed curves indicate isochores representing the lowest and highest density covered by the tables, respectively. As we have discussed earlier, the nonideal-mixing model we have used for Ga-IV/Fe in (b) stabilizes it at lower pressures and destabilizes at higher pressures, which is why the two solid solutions nearly swap stability fields from (a) to (b). Although not shown here, we have also produced diagrams for the case of 95% Fe and 5% Ga, which look qualitatively similar to those presented in this figure.

We may use the EOS tables produced by PSO to perform calculations, just like we can with any EOS from the LEOS or SESAME libraries mentioned in the Introduction. As an example, Figure compares Hugoniot and isentrope predictions from the Fe/Ga mixture tables with those from the pure Fe table by Benedict et al.[58] The initial pressure and density in all cases are 240 GPa and 12.15 g/cm3, respectively, which is in the fluid phase and could represent the result of a shock-melting experiment. The Hugoniot curves and isentropes in the figure might therefore represent double-shock or shock–ramp experiments, respectively. Interestingly, the effect of introducing Ga solute particles to the Fe solvent appears to be much more pronounced on the temperature than it is on the pressure or the shock speed–particle speed relation. Part of this is due to the ideal-mixing entropy terms of the form −Rz ln z, but nonideal contributions to the entropy and the heat capacity also have a noticeable effect. This suggests that obtaining reliable temperature measurements along dynamic-compression paths, which is an ongoing and extremely challenging area of work,[62] could be important toward discerning the effect of impurities under relevant high-pressure conditions.

Figure 9

Comparison of Hugoniot curves and isentropes of pure Fe vs those of Fe/Ga mixtures at two different overall compositions (90% Fe and 95% Fe) and described by different mixing models: (a–c) pressure, temperature, and shock speed along the Hugoniot and (d) temperature along the isentrope. (We have not included pressure–density results for the isentrope in the figure because the differences are even smaller than they are for the Hugoniot.) All of the curves are initiated at an initial pressure and density of 240 GPa and 12.15 g/cm3, respectively. As explained in the text, the results for Fe/Ga mixtures are calculated from the multiphase EOS tables for these mixtures that are produced from PSO-generated temperature–pressure phase diagrams like the ones shown in Figure .

Conclusions

Critical Evaluation of Our Approach and Future Work

Now that we have demonstrated the reliability and utility of PSO with the examples in the previous section, we are in a position to suggest future work to address the current limitations of our approach and expand its capabilities. We focus on limitations that we believe can be overcome in the near future. One topic that we do not dwell on at length here is chemical-reaction equilibria, which broadly includes ionic equilibria (salt solutions), acid–base equilibria, and other redox reactions that are of relevance to planetary science in particular. While PSO can describe chemical-reaction equilibria (e.g., Bonilla-Petriciolet and Segovia-Hernández[51] have applied it to model equilibria in reacting fluid mixtures at low pressures), we anticipate that for high-pressure applications, this will be more of a long-term challenge than the other topics that we describe below. The main difficulty is that reactions tend to occur not in isolation, but in a complex network that can involve hundreds of species linked through several different reactions. Each reaction, whose mechanism will likely not be well characterized, introduces additional mass-conservation constraints and constraints on the chemical potentials due to the requirement for vanishing entropy production at equilibrium.[63] These constraints will have to be verified, just like how we have verified that our PSO-predicted phase boundaries satisfy the equality of chemical potentials in eq . Moreover, EOS models will need to be supplied for each species in order to do a proper calculation of the equilibrium state. This will be a challenge for ionic and possibly also for acid–base equilibria because the charged species involved in them are notoriously difficult to model.[64] All of these difficulties will be exacerbated at high pressures because the application of pressure can introduce additional reactions and species in the complex network and cause the species to become more closely packed (and thereby interact more strongly), so that, for example, simple electrolyte models like Debye–Hückel that work well for dilute mixtures are not applicable. Thermochemical codes such as Magpie[13] and Cheetah[12,65−67] are able to model high-pressure, chemical-reaction equilibria of certain select mixtures (see section ), but they necessarily have to adopt some simplifications in order to make the problem tractable.

One way to extend our algorithm, which is currently limited to solutions, is to enable it to handle compounds of the type discussed in section . Because these compounds have a fixed stoichiometry, each possible phase of the compound brings in an additional phase mole fraction—but not an additional set of phase compositions—whose value is to be adjusted during the minimization process. For example, suppose that Fe3Ga is the first compound[40] to appear on the Fe-rich side of the Fe/Ga phase diagram and that it can exist in different possible phases. Then, for a given T and P and for all Fe/Ga mixtures where the Fe overall mole fraction is greater than 0.75, we can generalize the Gibbs energy that is to be minimized via PSO from its form in eq towhere the double sum over i and j concerns the three solution phases (fluid, Ga-III/Fe, and Ga-IV/Fe) considered in the previous section, and the summation over k is performed over the phases of Fe3Ga. Here, χ is the mole fraction of phase k of Fe3Ga, and G(T,P) is the Gibbs energy of that phase obtained from the chosen multiphase EOS for Fe3Ga. Clearly, this EOS must be constructed beforehand if it is to be used in eq , which speaks to the issue we raised earlier about the challenging need to construct separate multiphase EOSs for each compound of interest. In order to account for the presence of the compound Fe3Ga, the mole-balance constraints in eq must be modified to

We may apply PSO to minimize the Gibbs energy in eq subject to these mole-balance constraints to obtain the equilibrium state corresponding to the given T, P, and . The equilibrium state indicates the mole fractions χ and compositions z of the three solution phases and the mole fractions of the different phases of Fe3Ga. If the PSO algorithm is working correctly, the equilibrium mole fraction of all Fe3Ga phases except for that of the most stable phase(s) at the given T and P should be virtually zero, much like how PSO is able to select out and eliminate all unstable solution phases, as illustrated in the previous section. The same procedure may be followed for less than 0.75 with other compounds (such as Fe6Ga5, Fe3Ga4, and FeGa3), as long as multiphase EOSs are available for these compounds, so that we can compute their Gibbs energy.

A limitation of our PSO-based approach as presently designed is the computational efficiency. This is true in general for global optimization methods applied to phase-equilibria computations. For each set , the Fe/Ga examples we have presented take a few seconds to converge on a laptop equipped with a 2.8 GHz Intel Core i7 processor. The exact amount of time varies due to the stochastic nature of the optimization process, and it also depends on the number of phases present in the equilibrium state; fewer iterations are generally required to achieve convergence in single-phase regions than in two-phase regions. Construction of the phase diagrams in Figures –6 typically takes several hours because each diagram involves performing a Gibbs-energy minimization for thousands of different sets of . While the speed could certainly be improved significantly along the lines suggested below, PSO is still a practical means for Gibbs-energy minimization. It is significantly slower, however, for Helmholtz-energy minimization. Each set for the Fe/Ga examples takes hundreds of seconds to complete, which is 2 orders of magnitude slower than an equivalent calculation for a set in Gibbs-energy minimization. Part of the computational inefficiency of our code may be attributed to the fact that we are running a development version (not a production version) that is written in Python, an interpreted language. If it were instead to be optimized and rewritten in a compiled language like C++, it is reasonable to expect a several-fold (perhaps even more than a 10-fold) increase in efficiency. A more significant speedup could be achieved through parallelization in which the particles of a swarm are divided up among different processors. We note that PSO is amenable to architectures involving graphics processing units because the particles move largely independently of each other, with the only communication necessary being the need to track yswarm,best, which is the most optimal position in the search space ever accessed by any of the particles. Once the necessary speedup has been achieved, we aim to perform Helmholtz-energy minimization to construct temperature–density phase diagrams (for a fixed overall ) from which we can readily obtain EOS tables. This is a more direct way of producing these tables than the Gibbs-based procedure followed in this study, which requires post-processing steps to convert temperature–pressure phase diagrams to the desired tabular format.

Improving the computational efficiency will also be a necessary step toward developing a PSO capability for entropy maximization. The goal of entropy maximization is to find the equilibrium state of a multicomponent mixture for a given set of E, ρ, and , where E is the total molar internal energy of the mixture. The swarm in this case symbolizes a microcanonical ensemble because each PSO particle represents a mixture that is fixed at the specified E, ρ, and . The equilibrium state reveals the compositions of the solution phases and the mole fractions of all phases (both solutions and compounds) and the temperature and pressure that are common to all phases. Adapting the Gibbs-energy minimization to Helmholtz-energy minimization requires only one additional step; the same is true for entropy maximization, where the additional step in this case involves a simultaneous determination of both temperature and pressure (instead of only pressure like in Helmholtz-energy minimization) for each PSO particle at every iteration. While conceptually simple, this additional step will again be numerically expensive, even more so than the pressure-only solution required in Helmholtz-energy minimization. One possible way to ease the computational burden is to apply the scheme we have developed for a related problem, where we map our chosen set of free-energy models for the mixture onto a more simple EOS in which pressure depends linearly on volume.[68] The functional form of this simple EOS is specifically designed, so that the T and P consistent with the given E and ρ can be found analytically. This mapping onto the analytically invertible EOS would have to be performed for each PSO particle at every iteration, but we expect that it will be much more computationally efficient to do this than to carry out a direct numerical solve for T and P. The motivation for developing an entropy-maximization capability comes from the needs of the continuum-scale, multiphysics codes for high-pressure applications that we have mentioned in the Introduction. These codes track and update E, ρ, and along the different mesh points, but they do not track T and P or the phase mole fractions and compositions. Although this latter set of information is not tracked, the codes do need to extract this information from the given EOS models because it is of interest (especially T and P) to the various constitutive models in these codes, like those for strength effects, kinetics, transport coefficients, and so forth. If the computational efficiency can be sufficiently improved, it may be possible to implement a PSO-based, entropy-maximization procedure directly into these codes, so that they can perform on-the-fly, phase-equilibria calculations of multicomponent mixtures. This would be a thermodynamically rigorous alternative to simplified, albeit more pragmatic, approaches that treat these mixtures as if the components were homogeneously mixed (i.e., neglects the fact that they can separate into multiple phases with different compositions).

As the first demonstration in the literature of PSO on high-pressure phase equilibria, we have chosen to work with a binary system (Fe/Ga) of up to three phases because it is relatively easy to explain, validate, and visualize the performance of our method on such a system. Our algorithm is nonetheless designed to be applicable more generally to mixtures with an arbitrary number of components c and (solution) phases . Part of our future work will be to test its performance on true multicomponent mixtures. In anticipation of this future testing, we have performed a computational-cost analysis of Gibbs-energy minimization of ternary and quaternary mixtures. We have mimicked the additional components in the cost analysis by creating fictitious secondary Fe and Ga components with different molar masses. Table indicates an approximately linear scaling in both c and . The table gives only a rough idea of the scaling because the actual run time will vary due to the stochastic nature of the method and the complexity of the Gibbs-energy surface. It should also be pointed out that the table paints an incomplete picture because we have kept the swarm size Np fixed in it, yet the free-energy surface will become more complicated (will increase in dimensionality) as c increases, implying that Np will likely need to increase with c to ensure that the swarm successfully converges to the true minimum. Hence, we have performed a second cost analysis in which Np is varied while c and are kept fixed, the results of which are presented in Table . It indicates that the run time also scales linearly with Np. In summary, our analysis has demonstrated that the PSO performance cost scales linearly with c, , and Np, which is only a fairly modest increase. This suggests that PSO would continue to be a suitable approach for more complicated mixtures than those investigated in this study, especially taking into consideration the fact that PSO is amenable to parallelization.

Table 1

Computational Cost of PSO Calculations as a Function of the Number of Components c and Maximum Possible Number of Phases in the Mixture of Interesta

c = 1	1.00	1.35
c = 2	1.79	2.72	3.66
c = 3	2.89	4.32	5.82
c = 4	3.63	5.42	7.50

Each entry indicates the average time it takes to perform a Gibbs-energy minimization for a given relative to the simplest case of . Each listed time is an average over 20 PSO runs carried out with a fixed number of Np = 60 particles in the swarm. The case is not examined because it is forbidden by the Gibbs phase rule.[35−39,95]

Table 2

Computational Cost of PSO Calculations as a Function of Np for a Binary (c = 2) Mixture that can Exist in up to Phasesa

Np = 50	Np = 100	Np = 200	Np = 400	Np = 800
1.00	2.06	3.82	7.56	15.20

Each entry indicates the average time (averaged over 20 runs) it takes to perform a Gibbs-energy minimization for a given relative to the fastest case of Np = 50 particles.

Comparison with Other Approaches

Here, we briefly review three other approaches for modeling the thermodynamics of multiphase, multicomponent mixtures: (1) flash calculations,[69−77] (2) the thermochemical code Cheetah,[12,65−67] and (3) the CALculation of PHAse Diagrams (CALPHAD) methodology.[78−88] All three are now mature capabilities that have benefited from decades of development. We cannot make an apples-to-apples comparison because they each have certain functions that the others are not designed to do, but by pointing out some of these differences, we hope to clarify how this work fits in the context of the others.

The objective of a flash calculation is the same as that of our PSO-based optimization: determine the equilibrium state of a multicomponent mixture for a specified set of conditions given an EOS for these mixtures. The name stems from the widespread need in the chemical and petroleum industries to perform phase separations like “flash distillations” of fluid mixtures, and as such, virtually all of the published studies on this topic feature examples with fluid mixtures described by EOSs specialized for the low-pressure conditions (typically well below 1 GPa) of industrial relevance. Nevertheless, the thermodynamic principles that underpin flash calculations are also applicable to solid solutions and to higher-pressure conditions. The conventional approach to a flash calculation is to divide it into two sequential steps: a stability analysis that determines the number of stable phases, followed by a phase-split calculation that determines the phase compositions. Most of the algorithms have been developed for Gibbs-energy minimization,[69−71] though more recent work has also been done on Helmholtz-energy minimization[72,73] and entropy maximization.[74−77] The main advantage that these conventional flash methods offer over PSO is computational efficiency. This efficiency has allowed them to be implemented into the multiphysics simulation codes used by the chemical and petroleum industries. For Gibbs-energy minimization, a flash calculation tends to be several times faster than PSO. The disparity becomes even greater for Helmholtz-energy minimization and entropy maximization; the computational cost of a flash calculation for these more complicated problems increases by an order of magnitude or less, whereas it increases by at least 2 orders of magnitude for PSO. However, PSO offers distinct advantages that may make it a more attractive option under certain circumstances. First of all, it is far more suited to parallelization, which could potentially allow it to make up for a significant part of its extra computational requirements. It also does not perform the arbitrary separation into stability and phase-split steps. We have seen that going from Gibbs-energy minimization to Helmholtz-energy minimization to entropy maximization is conceptually a straightforward matter in PSO. In contrast, for flash calculations, these transformations often require a substantial reformulation of the algorithms. Perhaps most importantly, one does not need to supply PSO with a good initial guess for quantities like phase mole fractions and compositions (though the tradeoff for this robustness is the higher computational load). This is essential for high-pressure applications, as we have discussed in section . Flash calculations are the opposite in that they are very much reliant on a good initial guess for such quantities (this is especially true for Helmholtz-energy minimization and entropy maximization as explained in a recent study[77]), without which they are prone to failures like converging to an unphysical solution or not even converging to a solution at all. Correlations and recipes for initial guesses have been developed for a variety of fluid mixtures under low-pressure industrial conditions, but ultimately they are all based on trial-and-error and may not be applicable to higher pressures or to solid mixtures. As a result, it is not uncommon in a flash calculation to have to cycle through different initial guesses until finding one that works, which can significantly impact the robustness of these methods and add to the computational cost.

Cheetah[65,66] is widely used for modeling thermodynamics under detonation, which includes chemical-reaction equilibria. It does this by representing the interatomic interactions between the different species in a mixture with exponential-six potentials,[12,67] whose parameters are fit to experimental data. These interactions are translated to a free energy via the application of mixing rules and statistical–mechanical models, and the free energy is minimized to obtain the equilibrium state for a given set of conditions. Cheetah has been coupled with high-pressure continuum-scale codes, so that it can feed information directly to these codes. It can model systems ranging in complexity from monatomic materials to large organic molecules like high explosives and plastics. However, Cheetah cannot describe regimes where there is significant electronic ionization (recall from the Introduction that the LEOS and SESAME tables that we use for the pure components account for this through an electron-thermal term) and is limited in its applicable temperature and pressure range as a result. It is also not intended to model dense metals like pure Fe, pure Ga, and Fe/Ga mixtures because the exponential-six potentials are not expected to give accurate predictions for such materials over most conditions of interest.

CALPHAD[78−83] is the most well-known methodology for multiphase, multicomponent mixtures in materials science. It follows the standard prescription from chemical thermodynamics[35−39] outlined in section , in which the Gibbs energy of each solution phase is expressed as a linear combination of the Gibbs energies of the pure components that make up the phase plus some additional mixing terms. CALPHAD applies an optimization procedure (which could be a Newton-based method or a regression approach or a Bayesian technique[84]) to fit the parameters of the mixing models, and potentially also those of the pure-component Gibbs energies, to reproduce experimental data on the phase diagrams. Unlike Cheetah, it is designed to handle dense metallic systems, and in fact, producing phase diagrams for metallurgical applications is probably the most common use of CALPHAD. However, the vast majority of CALPHAD studies focus only on ambient pressure, where experimental phase-diagram data are far more readily available than at higher pressures. Some relatively recent studies have attempted to extend CALPHAD to higher pressures,[85−88] in which most of the pressure dependence is introduced by adding a cold energy and an ion-thermal term (an Einstein model with a pressure-dependent Einstein temperature), but again no electron-thermal term, to the ambient pure-component Gibbs energies. These high-pressure extensions are tied to the legacy of the earlier, ambient-pressure database through the use of an ad hoc interpolation function, and this constraint to the ambient-pressure models can in some cases make it difficult to obtain agreement with experimental data at higher pressures, such as those pertaining to melt curves.[87] While CALPHAD could be a good option for high-pressure applications, it is not widely used within high-pressure science at the moment because the details regarding the implementation of the underlying models and algorithms are obscured behind propietary software packages (e.g., Thermo-Calc, FactSage). Thus in practice, it requires an experienced user of these packages to make the modifications necessary to adapt the underlying databases to high pressures. Furthermore, to the best of our knowledge, the high-pressure CALPHAD-based models have remained as standalone objects in the sense that there has been no attempt to construct EOS tables from them or to otherwise connect CALPHAD to the continuum-scale codes mentioned in the Introduction.

Summary

This study is motivated by a number of practical needs for thermodynamic modeling of multicomponent mixtures at high pressures. One is the desire to have a framework for these mixtures that directly uses existing EOS models specifically designed for high-pressure conditions, with prominent examples being the EOSs that make up the LEOS and SESAME libraries. Section and the Appendix outline such a formulation in which mixing terms are added onto the pure-component properties that are obtained from the existing high-pressure EOSs. Our study is also motivated by the fact that multicomponent systems, like their single-component counterparts, can exist in a multiphase state when in equilibrium. A complete determination of the equilibrium state therefore requires finding the relative amounts of the different phases and the composition of each phase. This may be achieved through a free-energy minimization procedure. Such a capability would also be needed to generate phase diagrams because a phase diagram is merely a visual depiction of the equilibrium states of a mixture across a range of conditions. A crucial point that we reiterate here is that unlike the situation at ambient pressure, available data on phase diagrams at high pressures tend to be extremely limited. A free-energy minimization method that is suitable for high-pressure conditions must therefore be able to generate phase diagrams with little information (minimal input from the user regarding initial guesses, actual number of phases in equilibrium, etc.) on how the phase boundaries are connected together.

Our PSO-based minimization scheme described in section satisfies this important criterion. Another key reason why we have chosen PSO is because it is a simple method that does not require expertise in or access to any specialized software. As long as it is given free-energy models for each of the possible phases that might be present in the equilibrium state, it will iteratively adjust the relative amounts of these phases and the compositions of the individual solution phases to find the minimum on the Gibbs-energy surface provided by the models. In doing so, it will reveal the equilibrium state for a given temperature T, pressure P, and overall composition . If it turns out, for example, that a particular phase is not present in the equilibrium state at the specified , PSO will determine its mole fraction to be very close to zero, typically 10–9 or less. Section tests and demonstrates the reliability of our method with different example mixing models by generating temperature–composition phase diagrams. It is able to capture the fact that the particular nonideal-mixing models we have used for demonstration purposes in this study result in nontrivial features like eutectic and peritectic points. Moreover, we have applied PSO to produce EOS tables for mixtures that rigorously account for changes in phase boundaries due to mixing. As far as we know, this is the first demonstration of such a capability.

In addition to its reliable and readily accessible nature, PSO is also highly adaptable. We have discussed how our method, which is currently limited to mixtures involving fluid or solid (e.g., alloy) solutions, may be extended to also allow for stoichiometric compounds. This is a topic that we intend to work on in the near future. With minor modifications, the PSO algorithm may also be adapted to minimize the Helmholtz energy or maximize the entropy instead of minimizing the Gibbs energy. In fact, we have done this already for Helmholtz-energy minimization and have used it to cross-check results of the Gibbs-energy minimization. A limitation, however, is that our current Helmholtz-energy minimization procedure is too computationally inefficient to be used for its original intended purpose of constructing EOS tables in a direct manner through the generation of temperature–density phase diagrams. We have given suggestions on how to improve the efficiency in the previous section, though there are likely to be other avenues for speedup as well. We note that PSO is amenable to parallelization, and our tests show that it scales well (linearly) with respect to the number of components, phases, and particles in the swarm. With sufficient improvements in the efficiency, we hope to use PSO to produce multiphase EOS tables for mixtures via Helmholtz-energy minimization. It may also be possible to implement a PSO-based entropy-maximization procedure directly into continuum-scale codes for high-pressure applications to enable them to perform on-the-fly, phase-equilibria calculations of mixtures in a thermodynamically rigorous fashion.