QH-POCC: taming tiling entropy in thermal expansion calculations of disordered materials

While measuring CTEs has become an increasingly integral part of experimental research of MCAs and MCCs [12, 13, 14, 15], the substitutional disorder presents a formidable challenge for computational investigations of these materials, resulting in only few computed CTE values with scarce experimental validation [16, 17]. To date, calculations of thermal expansion in systems with substitutional disorder have employed the Coherent Potential Approximation (CPA) [18, 19], Special Quasirandom Structures (SQS) [20], or the cluster expansion (CE) method [21, 22]. The CPA allows for modeling alloys of any stoichiometry using relatively small unit cells. However, it cannot reliably calculate atomic forces [23] and is thus unable to employ conventional lattice dynamics methods. SQS models disordered systems at a fixed composition by finding the most random (infinite temperature) representation of the structure. As a consequence, it neglects finite-temperature short-range ordering effects. CE calculates thermodynamic properties through a series expansion of the free energy using a set of ordered structures. While the obtained free energy is exact, the combinatorial explosion of the number of cluster expansion parameters makes the study of complex MCAs and MCCs impractical.

The Partial OCCupation (POCC) method was developed to include finite-temperature effects while alleviating the computational cost of CE [24], thus allowing practical studies of multi-component systems. POCC models the disordered material through a series of small tiles — ordered representative supercells with the same stoichiometry as the parent structure — which are also capable of calculating vibrational properties [25]. This particular expansion allows POCC to factorize the tiling entropy associated with the latent heat of a first order phase transition (from chemical order and microstrucural disorder to chemical disorder and no microstrutures), and to generate a partial partition function out of the global one. The approach, taken at finite temperatures with the quasi-harmonic approximation (QHA), enables the calculation of thermomechanical properties. The extended method, called QH-POCC, is tested on AuCu₃ and CdMg₃ for the ordered and disordered fields. Both thermal expansion and isobaric heat capacity are in good agreement with experiments, demonstrating the potential of QH-POCC to accurately model thermomechanical properties in the QHA at a reasonable computational cost. We chose the latter compound for its capabilities in forming long-period superstructures at low temperature [26, 27, 28], so that the effect of the diverging characteristic lengths of the ground-state can be qualitatively compared against the finite size of the POCC tiles expansion.

Thermal Expansion of Disordered Systems. The volumetric CTE $\beta$ is defined as the logarithmic derivative of the temperature-dependent equilibrium volume:

where $V_{\textnormal{eq}}(T)$ is the equilibrium volume at temperature $T$, i.e., the volume $V$ for which the temperature- and volume-dependent free energy $F$ is at a minimum:

$F(V,T)$ consists of the following contributions when neglecting anharmonic effects:

$E_{0}$ is the energy of the structure at 0 K. $F_{\textnormal{elec}}$ is the electronic contribution to the free energy and can be calculated as:

Here, $g_{\textnormal{e}}\left(\epsilon,V\right)$ is the electronic density of states (DOS), $f_{\textnormal{FD}}=f_{\textnormal{FD}}\left(\epsilon,T\right)$ is the Fermi-Dirac distribution, $E_{\textnormal{\tiny F}}$ is the Fermi energy, and $s_{T}(\epsilon)=f_{\textnormal{FD}}\log(f_{\textnormal{FD}})+(1-f_{\textnormal{FD}})\log(1-f_{\textnormal{FD}})$ is the electronic entropy of $\epsilon$.

The last contribution, $F_{\textnormal{vib}}$, is the vibrational free energy. It can be calculated via the phonon DOS, $g_{\textnormal{ph}}\left(\omega(V)\right)$, as:

where $\omega$ is the phonon frequency and $\hbar$ is the reduced Planck constant.

In the QHA, $V_{\textnormal{eq}}(T)$ is obtained by calculating $F(V,T)$ over a set of volumes and fitting the volume-dependent free energy to an equation of state (EOS) for each temperature. Calculating the volume-dependency of $F$ is thus the central problem for the calculation of thermomechanical properties at finite temperatures.

Computational modeling of MCAs and MCCs — and therefore accurate calculation of free energies — is challenging due to the substitutional disorder in these compounds. AFLOW-POCC describes disordered materials through an ensemble of small Hermite-normal-form (HNF) ordered representatives [29, 30, 31], also called tiles. This method reproduces the correct composition, symmetry, and properties that depend on the radial distribution function of the system but needs to be integrated over all possible tiling configurations if the global partition function $Z_{\rm global}$ and free energy $F_{\rm global}$ are needed. For an $NVT$ ensemble, $Z_{\rm global}$ can be calculated as:

where $k_{\textnormal{\tiny B}}$ is the Boltzmann constant and $V_{i}$ is the volume of the $i^{\textnormal{th}}$ tile. The sum is constrained by the organization of the tiles covering the whole normalized volume $\sum\left\{V_{i}\right\}/N=V$ with the global number of tiles $N$. $Z_{\rm global}\left(V,T\right)$ can be used to obtain $F_{\rm global}$ through:

The calculation of $Z_{\rm global}\left(V,T\right)$ is highly non-trivial. The number of possible configurations $\Omega_{\rm tiling}$ — having the correct volume where all the pieces fit together — increases with the cut off of the maximum POCC tile-size. Thermodynamically, this number also leads to the tiling entropy $S_{\rm tiling}\equiv k_{\textnormal{\tiny B}}\log\Omega_{\rm tiling}$. The properties describing the “organization”, $\left\{\Omega_{\rm tiling},S_{\rm tiling}\right\}$, can be calculated combinatorially or through Monte Carlo modeling like in percolation or self-avoiding random-walks theories, but incurs great computational cost. Here, these quantities can be neglected because of the following considerations:
i. $\Omega_{\rm tiling}$ varies slowly unless near a first-order phase transition. The number of possible configurations depends on the identity and distribution of the individual tiles. Therefore, $\Omega_{\rm tiling}$ is nearly constant below and above the transition temperature, albeit not the same in general. At the phase transition, on the other hand, the Boltzmann population vector suddenly rotates in the probability hyperspace (see Fig. 2e in Ref. [32]), leading to completely different tiling-distribution identities and thus a different $\Omega_{\rm tiling}$. Similar considerations hold for the tiling entropy: $S_{\rm tiling}(T)$ is a slowly varying function of $T$ unless a phase transition is crossed.
ii. Derivatives $\partial[\cdot]/\partial T$ operating on functionals containing $S_{\rm tiling}(T)$ may ignore tiling contributions. As long as $\partial[\cdot]/\partial T$ is not operated at the phase transition (or nearby, as precursors start appearing [32]), $S_{\rm tiling}(T)$ can be neglected with great computational savings. Temperature-derivative properties, such as specific heat and thermal expansion, can be estimated below and above the transition temperature without direct knowledge of the organization of the tiles.
iii. Once the tiling organization is settled, phase stability requires minimization of free energies. This, in turn, requires minimization with respect the total tiling organization along the aforementioned volume constraints, $\sum\left\{V_{i}\right\}/N=V$ (Lie derivative). In general, $\min\left[\sum\left\{\cdot\right\}\right]\neq\sum\min\left[\left\{\cdot\right\}\right]$ — yet, the POCC model allows to solve this conundrum: in canonical POCC-tiles, both species concentration and superlattice are conserved. Only atoms are swapped generating symmetrically inequivalent decorations [24]. As such, the volumes $V_{i}$ and minimum energy $E_{i}$ (or $H_{i}$ in an applied stress field) for each $i$-tile will be quite close, $V_{i}\sim V$, so that the global EOS can be summarized as:

Equation (10) coupled with Equations (8) and (9) allow the definition of a partial partition function $Z_{\rm{\begin{subarray}{c}{POCC}\end{subarray}}}(V,T)$ and free energy $F_{\rm{\begin{subarray}{c}{POCC}\end{subarray}}}(V,T)$:

Even though $Z_{\rm{\begin{subarray}{c}{POCC}\end{subarray}}}$ and $F_{\rm{\begin{subarray}{c}{POCC}\end{subarray}}}$ neglect tiling entropy, they are capable of reproducing $T$-derivative properties outside regions of phase transition, without the cumbersome calculation of the global partition function. Equation (12) can then be used to construct the EOS of of the disordered material, leading to the global coefficient of thermal expansion via Equation (1). This scheme is illustrated in Fig. 1.

From the EOS-fitted curve, other thermomechanical properties can be derived. The isobaric ($C_{\textnormal{\tiny P}}$) and isochoric ($C_{\textnormal{\tiny V}}$) heat capacities can be calculated as:

where $B$ is the bulk modulus:

These quantities can be used to determine the effective Grüneisen parameter $\bar{\gamma}$:

Calculations’ details. AFLOW [40] leverages the Vienna ab initio simulation package (VASP) with all calculation parameters following the AFLOW standard [41]. Exchange and correlation were treated with the projector augmented wave (PAW) method [42] in either the local density approximation (LDA) [43] or generalized gradient approximation (GGA) proposed by Perdew, Burke, and Ernzerhof (PBE) [44]. The cutoff energies are chosen to be 1.4 times the recommended maximum cutoffs (ENMAX) of all pseudopotentials as set by VASP.

Disordered compounds were modeled using the AFLOW-POCC module with a supercell size of 4, resulting in 7 and 29 unique tiles, respectively. Phonon calculations were performed using the AFLOW-APL module [25]. Supercell sizes and $\mathbf{k}$-point mesh dimensions were selected such that free energies were converged to within 1 meV/atom for representative tiles of the disordered materials. For ordered representatives of CdMg₃, the supercells used to calculate the phonon properties were constructed to have at least 64 atoms, while for the ordered representatives of AuCu₃, the sizes were chosen as specified in Table 1. Forces within these supercells were calculated with a $\mathbf{k}$-points mesh having at least 6,900 and 3,000 $\mathbf{k}$-points per reciprocal atom (KPPRA) for AuCu₃ and CdMg₃, respectively [41].

For the EOS fit, vibrational properties were calculated at various compressed and expanded volumes, typically spanning -12% to 12% of the equilibrium volume at a 6% increment. For AuCu₃, adjustments needed to be made to accommodate for the presence of imaginary frequencies (see Table 1). Ordered representative #3 of AuCu₃ has unstable modes with the LDA functional at 3% volume compression, but these modes were easily omitted in the calculation of the thermomechanical properties because their contributions to the phonon DOS were negligible. Ordered representative #4 of AuCu₃ was found to be dynamically unstable and was discarded entirely from the ensemble of structures. The calculated free energies were fitted using the Stabilized Jellium EOS [45, 46]:

The parameters $f_{i}$ were determined using a polynomial fit.

Fig. 2 shows the calculated thermomechanical properties for these materials in their ordered and disordered phases using the Local Density Approximation (LDA) and the Perdew-Burke-Ernzerhof (PBE) functional. For both ordered and disordered AuCu₃, $\beta$ increases with increasing temperature until they plateau above 200 K (Fig. 2a). The CTE of the disordered structure is larger than of the ordered phase for all temperatures, which is consistent with experiment. PBE overestimates experimental $\beta$ values for both states. The disordered phase even appears to diverge at higher temperatures, which is expected for QHA, as anharmonic contributions become increasingly dominant. LDA shows excellent agreement with experimental values, suggesting that the choice of the functional is the reason for the bad agreement for PBE.

The CTE of CdMg₃ is shown in Figure 2b. The curves of the ordered phases show similar behavior as those for AuCu₃, except for $\beta$ for LDA, which converges towards the value of the disordered phase. Good experimental values for the ordered phase are unavailable in literature, but extrapolating from the tail found in the experimental values, it is anticipated to see good agreement. Above the experimental order-disorder transition temperature, the disordered phase shows good agreement for PBE whereas LDA underestimates $\beta$. This shows that the functional choice is critical when calculating thermal expansion of these materials. QH-POCC also shows a striking difference in the disordered phase compared to AuCu₃: $\beta$ of disordered CdMg₃ exhibits a peak at 50 K and 70 K for PBE and LDA, respectively.

The isobaric heat capacities shown in Fig. 2c for AuCu₃ and in Fig. 2d for CdMg₃ are less sensitive to the chosen functional. The agreement with experimental values for the ordered structures are excellent for both functionals. $C_{\textnormal{\tiny P}}$ for the disordered phases are underestimated, but still in good agreement with experiments. A peak can be observed in QH-POCC for both LDA and PBE in disordered CdMg₃, with the peak for LDA being at a higher temperature.

For the bulk modulus, experimental data is only available for AuCu₃ [38]. As Fig. 2e shows, PBE underestimates $B$ whereas LDA strongly overestimates it. $B$ is very sensitive to the chosen functional since it depends on the volume of the unit cell. This is not a shortcoming of QH-POCC, but of QHA, as this behavior was also found for simple ordered materials such as MgO and CaO [52]. To further substantiate that QH-POCC is consistent with other methods, $B$ was calculated at 0 K using the AFLOW Automatic Elastic Library (AEL) [53] — Fig. 2e demonstrates that AEL and QH-POCC are consistent with each other. The values for the disordered phase is smaller than for the ordered phase at all temperatures, which is consistent with experiments. However, $B$ drops much faster with temperature than experiments suggest, which has also been observed in ordered materials [52].

As for AuCu₃, CdMg₃ has a smaller bulk modulus in the disordered than in the ordered state, with LDA resulting in larger moduli than PBE. The temperature dependence of disordered CdMg₃ differs from AuCu₃ by showing an inflection, the end point of which coincides with the temperature of the peak in $\beta$ and $C_{\textnormal{\tiny P}}$ for both functionals.

Finally, Figs. 2g and h show the effective Grüneisen parameters for both materials. In general, $\bar{\gamma}$ is larger for the disordered material than the ordered material, particularly at higher temperatures. An exception is LDA and CdMg₃, where $\bar{\gamma}$ for the ordered phase approaches the disordered phase due to the divergence of $\beta$ at higher temperatures. Higher values of $\bar{\gamma}$ correlate with lower thermal conductivity, which is expected for materials with substitutional disorder [3]. The divergences found at low temperatures are the result of the calculation method: as Equation 16 shows, $\bar{\gamma}$ depends on the quotient of $\beta$ and $C_{\textnormal{\tiny V}}$, both of which are very small at low temperatures, causing numerical instabilities. This relationship also causes the peak for CdMg₃ as $\beta$ peaks stronger than $C_{\textnormal{\tiny P}}$.

In experiments, peaks in thermal expansion and isobaric heat capacity typically indicate phase transitions. Here, peaks are not found in QHA for ordered phases and in disordered CdMg₃ regions. From the solid solution and upon reducing temperature, CdMg₃ transitions from a disordered hexagonal close-packed structure to the ordered D0₁₉ [26, 51]. With further temperature reduction, D0₁₉ often develops antiphase boundaries and complex stacking faults [39], associated with forming long-period superstructures, LPS. Notably, several Cd and Mg binary alloys show similar features [26, 27, 28]. The divergence of the characteristic length scale in the LPS can not be reproduced by the tiles having a finite size, making the POCC expansion incomplete. Therefore, the constrained LTVC population vector [32] will rotate and promote modifications of D0₁₉ containing antiphase boundaries and stacking faults having increased volume without ever reaching full LPS description. While such computational limitation is expected to appear at low temperatures and in all systems having ordered ground-states with characteristic length scales too large for the POCC-tiles expansion, our algorithm has been devised to study disordered systems at high temperatures where only short-range needs to be described [11, 54]. As such, the limitation does not pose any concern.

Experiments suggest that a peak should appear in $\beta$ and $C_{\textnormal{\tiny P}}$ for AuCu₃ as well, but none are observable in Figs. 2a and c. The reason for this discrepancy could be the sample size of the canonical ensemble: POCC creates only 7 ordered representatives for AuCu₃, but 29 for CdMg₃. To test this hypothesis, the size of the HNF matrix was increased to 8, which produced 49 ordered representatives. As Figs. 3a and b show, increasing the maximum POCC tile size results in a peak in both thermal expansion and heat capacity, confirming that small sampling suppressed the appearance of the phase transition. However, no inflection point can be observed in the bulk modulus as was observed in CdMg₃ (see Fig. 3c). This is not surprising because the values for the ordered and disordered phase are similar at low temperatures in both QH-POCC and experiments, making any inflection imperceptible. Finally, the effective Grüneisen parameter, shown in Fig. 3d, shows a deep dip at lower temperatures when the supercell size is increased due to the larger peak height in $C_{\textnormal{\tiny P}}$, which is in the denominator of $\bar{\gamma}$.

As with CdMg₃, the calculated transition temperature of AuCu₃ is far below the experimental value. While transitioning from an ordered to a disordered phase, both materials are in a partially ordered state spanning at least 100 K. The requirement of POCC that the stoichiometry of the tiles correspond to the disordered structure parametrizes this transition region as a line. To improve the accuracy of the transition temperature, a more complete picture of the tile distribution is needed so that the variation of $S_{\rm tiling}$ can be captured. This could be achieved, for example, by including structures with compositions different from the parent structure. A similar approach is employed in the Lederer-Toher-Vecchio-Curtarolo (LTVC) method, which can accurately predict the transition temperature, but requires more input data [32].

Outside the transition region, however, all properties are well converged with respect to the POCC supercell size. So, while $\Omega_{\rm tiling}$ and $S_{\rm tiling}$ may change with increased maximum tile size, their temperature derivatives do not. The smallest supercell is thus sufficient to calculate thermophysical properties of the disordered material.