Here we briefly review some elementary details of alloy thermodynamics and statistical mechanics. You can also find much of this information on the arXiv: arXiv:2505.05393.

Alloy Thermodynamics

In a substitutional alloy with fixed underlying lattice, a particular arrangement of atoms can be specified by a discrete set of site occupation numbers, \({ \xi_{i\gamma} }\), where \(\xi_{i \gamma} = 1\) if site \(i\) is occupied by an atom of species \(\gamma\) and \(\xi_{i \gamma} = 0\) otherwise. The lattice index \(i\) takes values in the range 1 to \(N\), where \(N\) is the total number of lattice sites in the system, while the species index \(\gamma\) takes values in the range 1 to \(s\), where \(s\) is the number of chemical species (elements) present in the alloy composition. The physical constraint that each lattice site is occupied by one (and only one) atom is expressed as \(\sum_{\gamma} \xi_{i\gamma} = 1\), while the total concentration of a given chemical species is given by \(c_\gamma = \frac{1}{N} \sum_i \xi_{i\gamma}\). (Naturally, vacancies can be treated in this framework by considering them as an additional chemical species present at a very low concentration.) If we consider an ensemble of alloy configurations, we can then define the site-wise concentrations, \(\{ c_{i\gamma}\}\), as the average value of the site occupation numbers across the ensemble, \(c_{i\gamma} = \langle \xi_{i\gamma} \rangle\), where \(\langle \cdot \rangle\) denotes the average taken over the ensemble. (The precise meaning of ‘ensemble’ will be defined presently.)

In the context of substitutional alloys, the word phase refers to a physically homogeneous area of the material with uniform chemical composition and physical properties. At the atomic level, this means that, on average, the site-wise concentrations in a particular spatial region are homogeneous and describe some repeating motif of atoms arranged on a lattice. In general, a particular single-phase region of an alloy can be either a ‘solid solution’, in which all lattice sites have equal partial lattice site occupancies and elements mix randomly and homogeneously, or an ordered intermetallic phase (sometimes referred to as an ‘intermetallic compound’), where a crystallographically ordered, repeating motif of elements can be identified. Naturally, an intermetallic phase can also admit substitutional disorder on one or more of the atomic sites in the repeating motif. In thermal equilibrium, it is possible for more than one phase to be present in an alloy, with the maximum permitted number of phases present determined, in general, by the Gibbs phase rule. An alloy which decomposes into multiple coexisting phases in thermal equilibrium is said to undergo ‘phase segregation’ or ‘phase decomposition’. (Sometimes this process can also be referred to as one or more phases ‘precipitating’ out of the solid solution.)

For a system with fixed underlying lattice, in the high temperature limit, entropy dictates that atoms should have no lattice site preference. However, as the temperature is lowered, this site symmetry will eventually be broken and either a partial or full site ordering (or decomposition) will become established. The relevant equilibrium phase(s) of an alloy are, in principle, fully determined once the pressure/volume, temperature, and alloy composition (i.e. set of alloying elements and their concentrations) have been specified. Freezing the underlying crystal lattice removes the degree of freedom provided by the pressure/volume, and it remains to inspect the probabilities of a given arrangement of atoms occurring at the specified temperature.

In thermal equilibrium, the relevant probability distribution for determining the likelihood of a given configuration occurring is the Boltzmann distribution. The partition function, \(Z\), is written as

\[ Z = \sum_{ \{ \xi_{i \gamma} \} } e^{-\beta E( \{ \xi_{i \gamma} \} )}, \]

where \(E( \{ \xi_{i \gamma} \} )\) is the energy associated with a configuration, \(\{ \xi_{i \gamma} \}\), and the sum is taken over all possible atomic configurations. The symbol \(\beta\) is defined via \(\beta=1/k_B T\), where \(T\) is temperature, and \(k_B\) is the usual Boltzmann constant. The probability of a given arrangement of atoms occurring is then

\[ P(\{ \xi_{i\gamma} \}) = \frac{e^{-\beta E( \{ \xi_{i \gamma} \} )}}{Z}. \]

This probability distribution defines an ensemble of configurations distributed according to that probability. If the distribution is known, it is possible to recover various thermodynamic quantities of interest, such as the free energy and specific heat.

However, for any physically realistic form of the alloy internal energy, \(E( \{ \xi_{i \gamma} \} )\), and for all but the smallest of simulation supercells, direct evaluation of the partition function is computationally intractable due to the huge number of configurations which must be considered. This is a particular problem in the context of alloys containing multiple elements, such as high-entropy alloys—those alloys containing four or more elements alloyed in near-equal ratios—as the size of the configuration space grows combinatorially both with the total number of atoms in the simulation cell, as well as with the number of elements present in a given composition. Consequently, it is necessary both to find means by which to evaluate the energy associated with a given arrangement of atoms which are accurate and computationally efficient, as well as to use sampling algorithms to reliably estimate thermodynamic quantities and determine equilibrium phases as a function of temperature without evaluating the partition function in a brute-force manner. Evaluation of the internal energy of the alloy will be discussed now, while specific sampling algorithms (those which are implemented in BraWl) will be discussed later in this article.

Evaluation of the alloy internal energy

In the context of simulations performed on alloy supercells, the energy of a given configuration can be evaluated in a variety of ways. These include first-principles electronic structure calculations using density functional theory (DFT); interatomic potentials, including machine-learned interatomic potentials (MLIPs); and lattice-based models such as cluster expansions (CEs). All of these methods have their advantages and disadvantages. DFT calculations on alloy supercells, while highly accurate, are computationally expensive, rendering this option inviable when a large number of alloy energy evaluations are required for sampling. Calculations using interatomic potentials, MLIPs and CEs are significantly cheaper than direct DFT calculations, but they still frequently require a large DFT training dataset. Additionally, the number of fitting parameters required for these models grows significantly once a multicomponent alloy is considered, making them challenging to train and leading to concerns regarding their accuracy. Finally, it should be stressed that the computational cost of these models often remains prohibitive when a large number of energy evaluations is required.

An alternative, physically intuitive model for the internal energy of an alloy is the Bragg-Williams model, which assumes that the internal energy of an alloy takes a simple, pairwise form. The Bragg-Williams {\color{blue} model} Hamiltonian has the form

\[ H(\{\xi_{i\gamma}\}) = \frac{1}{2}\sum_{i \gamma; j\gamma'} V_{i\gamma; j\gamma'} \xi_{i \gamma} \xi_{j \gamma'}, \]

where \(V_{i\gamma; j\gamma'}\) denotes the effective pair interaction (EPI) between an atom of chemical species \(\gamma\) on lattice site \(i\) and an atom of chemical species \(\gamma'\) on lattice site \(j\). (The factor of \(\frac{1}{2}\) eliminates double-counting in the summation.) For a system of finite size, it is assumed that periodic boundary conditions are applied in all coordinate directions. We note that the form of this Hamiltonian is very similar to that of the Lenz-Ising model, an elementary model in magnetism.

Generally, the assumption is made that that the EPIs are spatially homogeneous and isotropic, and Eq.~eq:b-w1} is therefore rewritten as

\[ H(\{\xi_{i\gamma}\}) = \frac{1}{2}\sum_{i \gamma} \xi_{i \gamma} ( \sum_{n} \sum_{j \in n(i)} \sum_{\gamma'} V^{(n)}_{\gamma \gamma'} \xi_{j \gamma'} ), \]

where the sum over \(i\) remains a sum over lattice sites, but the sum over \(n\) denotes a sum over the coordination shells (nearest-neighbours, next-nearest-neighbours, etc.) of the lattice. The notation \(n(i)\) is then used to denote the set of lattice sites which are \(n\)th nearest-neighbours to site \(i\). Then \(V^{(n)}_{\gamma \gamma'}\) denotes the effective pair interaction between chemical species \(\gamma\) and \(\gamma'\) on coordination shell \(n\). It is reasonable to assume that, for most alloys, the strength of EPIs will tail off quickly with decreasing distance, and the sum over \(n\) can be taken over the first few coordination shells of the underlying lattice type being considered. (This is, of course, equivalent to imposing some radial ‘cutoff’ on an interatomic potential.)

EPIs for the Bragg-Williams Hamiltonian can be obtained using a variety of methods, generally those based on density functional theory. Similarly to the CE method, it is naturally possible to fit EPIs for a given alloy composition to a set of DFT total energy evaluations on alloy supercells with success. However, most frequently, such EPIs are obtained using the Korringa–Kohn–Rostoker (KKR) formulation of density functional theory, where the coherent potential approximation (CPA) can be used to describe the average electronic structure and consequent internal energy of the disordered alloy. There are then a variety of suitable techniques available for assessing the energetic cost of applied, inhomogeneous chemical perturbations to the CPA reference medium which naturally lead to extraction of EPIs. Such techniques include both the generalised perturbation method (GPM), as well as techniques using the language of concentration waves to describe the atomic-scale chemical fluctuations. Approaches based on concentration waves have been derived for alloys both in the binary and multicomponent settings. Once the EPIs for a given alloy composition are obtained, the phase stability of a particular alloy can be examined using sampling techniques applied to the Bragg-Williams model. This is the purpose of BraWl as presented in this work.

Model limitations

It should be emphasised that there are two key limitations to the above discussion, and to the applicability of the Bragg–Williams model more generally.

The first limitation is the lack of consideration of entropic contributions beyond that made by the configurational entropy, such as vibrational, magnetic, and electronic entropies. Perhaps most important is the consideration of vibrational entropy, which is most relevant when a system transitions between two phases with different elastic properties. For example, if (upon cooling) a system undergoes a phase transition from a disordered phase to an ordered phase, the latter of which is elastically stiffer than the former, the difference in vibrational entropy between the two phases can thermodynamically stabilise the disordered phase and lower the disorder-order transition temperature. Such effects are not accounted for in the Bragg-Williams model directly. However, they can be included in some approximate way by modifying EPIs to account for these effects. More generally, the effect of other entropic contributions, such as magnetic and electronic contributions, should also be considered as appropriate to a given system of study.

The second limitation is the assumption of a fixed underlying crystal lattice. In an alloy where there are atomic size mismatches, there will often be local distortions to the underlying crystal lattice to accommodate the mismatches. (This is because configurations with such local lattice distortions represent true minima of the potential energy surface.) It is understood that these effects can affect predicted disorder-order transition temperatures in alloys.

Neither of the above effects (additional entropic contributions or local lattice distortions) are explicitly considered in the fixed-lattice Bragg-Williams model as implemented in BraWl, and users should therefore take due care when interpreting simulation results.