Zero Static Internal Stress Approximation (ZSISA)

• ZSISA is an approximation within the quasi-harmonic framework that decouples internal coordinate relaxations from external strain optimization to simplify free-energy minimization.
• The method first relaxes atomic positions at the static Born–Oppenheimer level and then computes phonons on these fixed geometries, significantly reducing computational costs.
• ZSISA accurately predicts equilibrium lattice parameters and elastic constants but may underperform when internal relaxations strongly influence thermal responses.

Searching arXiv for papers on ZSISA and related quasi-harmonic methods. The Zero Static Internal Stress Approximation (ZSISA) is an approximation within the quasi-harmonic approximation (QHA) for determining finite-temperature equilibrium structures and derived response functions while avoiding the full cost of minimizing the free energy with respect to both external strains and internal coordinates at every temperature and stress state. In its standard formulation, the total free energy is written as $$F(\xi,T)=U(\xi)+F_{\rm ph}(\xi,T)+F_{\rm el}(\xi,T)$$, or equivalently $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$, and ZSISA replaces the full free-energy stationarity conditions by a two-stage procedure: internal coordinates are first relaxed only at the static Born–Oppenheimer level for each externally strained geometry, and phonons are then computed on that zero-static-stress manifold; the final minimization is performed only over the external structural variables (Gong et al., 21 Jul 2025). The approximation has become important because full anisotropic QHA and full free-energy minimization (FFEM) can require several dozens or even hundreds of phonon calculations, whereas ZSISA and its volume-constrained variants reduce this burden substantially while often preserving high accuracy for equilibrium lattice parameters, elastic constants, and thermal expansion (Rostami et al., 12 Mar 2025).

1. Definition within quasi-harmonic free-energy minimization

In QHA, the equilibrium structure at temperature $T$ and, when relevant, under external stress or pressure, is obtained by minimizing the Helmholtz or Gibbs free energy with respect to all structural degrees of freedom. For hcp metals this may be expressed in terms of $\xi=(a,c/a)$, whereas in the more general crystallographic formulation one uses $\{C_\gamma\}$ to denote lattice lengths, angles, and internal atomic coordinates (Gong et al., 21 Jul 2025). The free energy contains a static electronic term and a vibrational term, with the phonon contribution written as

$$F_{\rm ph}(\xi,T)= \frac{1}{2N}\sum_{\mathbf q,\eta}\hbar\,\omega_{\eta}(\mathbf q,\xi)+ \frac{1}{N\beta}\sum_{\mathbf q,\eta} \ln\bigl[1-e^{-\beta\hbar\omega_{\eta}(\mathbf q,\xi)}\bigr], \qquad \beta=(k_B T)^{-1},$$

or equivalently as

$$F_{\rm vib}(C,T)= k_B T \sum_{q\nu}\ln\! \bigl[1-e^{-\hbar\omega_{q\nu}(C)/k_B T}\bigr] +\frac12\sum_{q\nu}\hbar\omega_{q\nu}(C).$$

The computational bottleneck is that phonon spectra must, in principle, be evaluated across a multidimensional structural grid (Rostami et al., 2024).

ZSISA reduces the dimensionality of this minimization. In one formulation, for each applied external strain state one relaxes internal coordinates so that the static Born–Oppenheimer forces vanish, computes phonons only on these partially relaxed structures, and then minimizes the total free energy with respect to the external strains while ignoring phonon-induced internal forces (Rostami et al., 2024). In the osmium study this is described as decoupling the minimization over lattice strains from that over internal ionic coordinates by assuming that the static internal stress vanishes at each strained geometry; the resulting approximate free energy is

$$F_{\rm ZSISA}(\varepsilon,T) = U\bigl(\xi(\varepsilon)\bigr) + F_{\rm ph}\bigl(\xi(\varepsilon),T\bigr) + F_{\rm el}\bigl(\xi(\varepsilon),T\bigr),$$

where $\xi(\varepsilon)$ lies on the zero-static-stress manifold determined at each strain (Gong et al., 21 Jul 2025).

The practical meaning of the approximation is that the vibrational free energy contributes to the optimization only through the external variables. The internal coordinates are not re-optimized with respect to $F_{\rm ph}$ or $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$0 at each temperature. This is the origin of both the efficiency and the limitations of ZSISA.

2. Variational structure and the volume-constrained form

The exact QHA minimization requires simultaneous stationarity with respect to external and internal variables. In the notation used for polar and anisotropic crystals, one writes the structural variables as $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$1 or, more generally, $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$2, and the exact conditions are

$$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$3

ZSISA replaces the second condition by the static constraint

$$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$4

or equivalently

$$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$5

for each internal degree of freedom at fixed external parameters (Rostami et al., 2024). One then minimizes only over the external strain coordinates.

A widely used specialization is the volume-constrained variant, denoted V-ZSISA or v-ZSISA. In the volumetric form, one regards the volume $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$6 as the sole external parameter, relaxes all other lattice-shape variables and internal coordinates at each fixed volume, and minimizes

$$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$7

through the zero-pressure condition

$$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$8

This restricts the sampling of phonons to a one-dimensional path in structure space (Rostami et al., 2024).

For hcp metals, the same idea is often stated in terms of the $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$9 “stress-pressure” curve in the $T$0 plane. One defines $T$1 by

$T$2

and then evaluates

$T$3

In this approximation, thermal expansion is assumed to occur only through the unit-cell volume, while the lattice shape follows the $T$4 stress-pressure curve (Gong et al., 21 Jul 2025).

This distinction is central. ZSISA retains anisotropic coupling among external lattice degrees of freedom, whereas v-ZSISA suppresses that coupling by collapsing the optimization onto a single volumetric coordinate. A common misconception is that these are interchangeable. The available studies show that they can be close for some hcp metals, yet materially different when anisotropy is itself the target quantity (Rostami et al., 12 Mar 2025).

3. Accuracy for hcp metals: osmium and beryllium

The most detailed benchmarks in the provided literature concern hcp osmium and hcp beryllium. In osmium, the QHA study reports negligible deviations between ZSISA and FFEM, and also states that V-ZSISA influences the results very little (Gong et al., 21 Jul 2025). Quantitatively, with the PBE functional, static relaxation changes $T$5 by only $T$6 kbar $T$7 and $T$8 by $T$9 kbar $\xi=(a,c/a)$0. A direct ZSISA-versus-FFEM comparison at $\xi=(a,c/a)$1 K gives deviations below $\xi=(a,c/a)$2 kbar in each isothermal elastic constant, corresponding to better than $\xi=(a,c/a)$3. The difference between V-ZSISA along the $\xi=(a,c/a)$4 stress-pressure curve and full-grid interpolation along the $\xi=(a,c/a)$5 kbar isobar is likewise reported as $\xi=(a,c/a)$6 kbar for the critical elastic constants. The study concludes that the combined ZSISA+V-ZSISA errors never exceed $\xi=(a,c/a)$7 in elastic constants or in equilibrium lattice parameters (Gong et al., 21 Jul 2025).

The beryllium study reaches a closely related conclusion while making the comparison to FFEM explicit over temperature. At $\xi=(a,c/a)$8 K, the differences between ZSISA and FFEM are reported as $\xi=(a,c/a)$9 kbar $\{C_\gamma\}$0 and $\{C_\gamma\}$1 kbar $\{C_\gamma\}$2. At $\{C_\gamma\}$3 K, the deviations become $\{C_\gamma\}$4 kbar $\{C_\gamma\}$5 and $\{C_\gamma\}$6 kbar $\{C_\gamma\}$7 (Gong et al., 2024). For the volume-constrained comparison in the quasi-static framework, the study finds that V-ZSISA almost fully captures the effect of thermal expansion on $\{C_\gamma\}$8 and $\{C_\gamma\}$9.

These results establish a characteristic regime in which ZSISA is reliable: hcp metals with relatively small internal relaxations under strain and modest coupling between those internal displacements and phonons. The beryllium analysis attributes the success of ZSISA to the fact that the only required internal relaxation under uniaxial or bi-axial strain is a small shift $$F_{\rm ph}(\xi,T)= \frac{1}{2N}\sum_{\mathbf q,\eta}\hbar\,\omega_{\eta}(\mathbf q,\xi)+ \frac{1}{N\beta}\sum_{\mathbf q,\eta} \ln\bigl[1-e^{-\beta\hbar\omega_{\eta}(\mathbf q,\xi)}\bigr], \qquad \beta=(k_B T)^{-1},$$0 of the two atoms along the $$F_{\rm ph}(\xi,T)= \frac{1}{2N}\sum_{\mathbf q,\eta}\hbar\,\omega_{\eta}(\mathbf q,\xi)+ \frac{1}{N\beta}\sum_{\mathbf q,\eta} \ln\bigl[1-e^{-\beta\hbar\omega_{\eta}(\mathbf q,\xi)}\bigr], \qquad \beta=(k_B T)^{-1},$$1-axis, and that the coupling constant entering $$F_{\rm ph}(\xi,T)= \frac{1}{2N}\sum_{\mathbf q,\eta}\hbar\,\omega_{\eta}(\mathbf q,\xi)+ \frac{1}{N\beta}\sum_{\mathbf q,\eta} \ln\bigl[1-e^{-\beta\hbar\omega_{\eta}(\mathbf q,\xi)}\bigr], \qquad \beta=(k_B T)^{-1},$$2 is modest (Gong et al., 2024). The osmium study expresses the same point in more general variational language: if the phonon-induced internal stress is much smaller than the static stress arising from unrelaxed ions, and if internal-relaxation displacements are small at each strain, the ZSISA error is formally of higher order and numerically typically below $$F_{\rm ph}(\xi,T)= \frac{1}{2N}\sum_{\mathbf q,\eta}\hbar\,\omega_{\eta}(\mathbf q,\xi)+ \frac{1}{N\beta}\sum_{\mathbf q,\eta} \ln\bigl[1-e^{-\beta\hbar\omega_{\eta}(\mathbf q,\xi)}\bigr], \qquad \beta=(k_B T)^{-1},$$3 for the elastic constants of hcp Os (Gong et al., 21 Jul 2025).

A second important benchmark concerns ZSISA versus the quasi-static approximation (QSA). Both hcp studies report that QHA-derived elastic constants agree much better with experiment than QSA results, indicating that ZSISA’s main competition is not FFEM but rather the omission of vibrational free-energy effects altogether (Gong et al., 21 Jul 2025).

4. Efficient anisotropic implementations and Taylor-expanded vibrational free energies

Although ZSISA is already cheaper than FFEM, its direct anisotropic implementation can remain expensive because one may still need phonons on a multidimensional strain grid. Two of the cited works address this by combining ZSISA with low-order Taylor expansions of the vibrational free energy (Rostami et al., 2024).

In the anisotropic formulation, the vibrational free energy is expanded around a reference strain $$F_{\rm ph}(\xi,T)= \frac{1}{2N}\sum_{\mathbf q,\eta}\hbar\,\omega_{\eta}(\mathbf q,\xi)+ \frac{1}{N\beta}\sum_{\mathbf q,\eta} \ln\bigl[1-e^{-\beta\hbar\omega_{\eta}(\mathbf q,\xi)}\bigr], \qquad \beta=(k_B T)^{-1},$$4 up to second order:

$$F_{\rm ph}(\xi,T)= \frac{1}{2N}\sum_{\mathbf q,\eta}\hbar\,\omega_{\eta}(\mathbf q,\xi)+ \frac{1}{N\beta}\sum_{\mathbf q,\eta} \ln\bigl[1-e^{-\beta\hbar\omega_{\eta}(\mathbf q,\xi)}\bigr], \qquad \beta=(k_B T)^{-1},$$5

with $$F_{\rm ph}(\xi,T)= \frac{1}{2N}\sum_{\mathbf q,\eta}\hbar\,\omega_{\eta}(\mathbf q,\xi)+ \frac{1}{N\beta}\sum_{\mathbf q,\eta} \ln\bigl[1-e^{-\beta\hbar\omega_{\eta}(\mathbf q,\xi)}\bigr], \qquad \beta=(k_B T)^{-1},$$6. The corresponding vibrational stress is then obtained from the first and second strain derivatives, and the ZSISA equilibrium condition is written as

$$F_{\rm ph}(\xi,T)= \frac{1}{2N}\sum_{\mathbf q,\eta}\hbar\,\omega_{\eta}(\mathbf q,\xi)+ \frac{1}{N\beta}\sum_{\mathbf q,\eta} \ln\bigl[1-e^{-\beta\hbar\omega_{\eta}(\mathbf q,\xi)}\bigr], \qquad \beta=(k_B T)^{-1},$$7

This converts the problem from a multidimensional free-energy surface minimization to the solution of a small set of linearized equations (Rostami et al., 12 Mar 2025).

The computational consequences are substantial. For the “ZSISA–E$$F_{\rm ph}(\xi,T)= \frac{1}{2N}\sum_{\mathbf q,\eta}\hbar\,\omega_{\eta}(\mathbf q,\xi)+ \frac{1}{N\beta}\sum_{\mathbf q,\eta} \ln\bigl[1-e^{-\beta\hbar\omega_{\eta}(\mathbf q,\xi)}\bigr], \qquad \beta=(k_B T)^{-1},$$8Vib2” implementation, the number of phonon band-structure calculations is

• $$F_{\rm ph}(\xi,T)= \frac{1}{2N}\sum_{\mathbf q,\eta}\hbar\,\omega_{\eta}(\mathbf q,\xi)+ \frac{1}{N\beta}\sum_{\mathbf q,\eta} \ln\bigl[1-e^{-\beta\hbar\omega_{\eta}(\mathbf q,\xi)}\bigr], \qquad \beta=(k_B T)^{-1},$$9 for cubic systems,
• $$F_{\rm vib}(C,T)= k_B T \sum_{q\nu}\ln\! \bigl[1-e^{-\hbar\omega_{q\nu}(C)/k_B T}\bigr] +\frac12\sum_{q\nu}\hbar\omega_{q\nu}(C).$$0 for hexagonal, trigonal, and tetragonal systems,
• $$F_{\rm vib}(C,T)= k_B T \sum_{q\nu}\ln\! \bigl[1-e^{-\hbar\omega_{q\nu}(C)/k_B T}\bigr] +\frac12\sum_{q\nu}\hbar\omega_{q\nu}(C).$$1 for orthorhombic systems,
• $$F_{\rm vib}(C,T)= k_B T \sum_{q\nu}\ln\! \bigl[1-e^{-\hbar\omega_{q\nu}(C)/k_B T}\bigr] +\frac12\sum_{q\nu}\hbar\omega_{q\nu}(C).$$2 for monoclinic systems, and
• $$F_{\rm vib}(C,T)= k_B T \sum_{q\nu}\ln\! \bigl[1-e^{-\hbar\omega_{q\nu}(C)/k_B T}\bigr] +\frac12\sum_{q\nu}\hbar\omega_{q\nu}(C).$$3 for triclinic systems (Rostami et al., 12 Mar 2025).

The volumetric study develops an analogous hierarchy for $$F_{\rm vib}(C,T)= k_B T \sum_{q\nu}\ln\! \bigl[1-e^{-\hbar\omega_{q\nu}(C)/k_B T}\bigr] +\frac12\sum_{q\nu}\hbar\omega_{q\nu}(C).$$4:

• E2Vib1: Born–Oppenheimer energy to second order, $$F_{\rm vib}(C,T)= k_B T \sum_{q\nu}\ln\! \bigl[1-e^{-\hbar\omega_{q\nu}(C)/k_B T}\bigr] +\frac12\sum_{q\nu}\hbar\omega_{q\nu}(C).$$5 to first order; two phonon spectra.
• E$$F_{\rm vib}(C,T)= k_B T \sum_{q\nu}\ln\! \bigl[1-e^{-\hbar\omega_{q\nu}(C)/k_B T}\bigr] +\frac12\sum_{q\nu}\hbar\omega_{q\nu}(C).$$6Vib1: exact Born–Oppenheimer energy, $$F_{\rm vib}(C,T)= k_B T \sum_{q\nu}\ln\! \bigl[1-e^{-\hbar\omega_{q\nu}(C)/k_B T}\bigr] +\frac12\sum_{q\nu}\hbar\omega_{q\nu}(C).$$7 to first order; two spectra.
• E$$F_{\rm vib}(C,T)= k_B T \sum_{q\nu}\ln\! \bigl[1-e^{-\hbar\omega_{q\nu}(C)/k_B T}\bigr] +\frac12\sum_{q\nu}\hbar\omega_{q\nu}(C).$$8Vib2: $$F_{\rm vib}(C,T)= k_B T \sum_{q\nu}\ln\! \bigl[1-e^{-\hbar\omega_{q\nu}(C)/k_B T}\bigr] +\frac12\sum_{q\nu}\hbar\omega_{q\nu}(C).$$9 to second order; three spectra.
• E$$F_{\rm ZSISA}(\varepsilon,T) = U\bigl(\xi(\varepsilon)\bigr) + F_{\rm ph}\bigl(\xi(\varepsilon),T\bigr) + F_{\rm el}\bigl(\xi(\varepsilon),T\bigr),$$0Vib4: $$F_{\rm ZSISA}(\varepsilon,T) = U\bigl(\xi(\varepsilon)\bigr) + F_{\rm ph}\bigl(\xi(\varepsilon),T\bigr) + F_{\rm el}\bigl(\xi(\varepsilon),T\bigr),$$1 to fourth order; five spectra (Rostami et al., 2024).

Across twelve benchmark materials, the volumetric study reports that three full phonon spectra, corresponding to quadratic order, are enough to determine the thermal expansion in reasonable agreement with the v-ZSISA-QHA method up to $$F_{\rm ZSISA}(\varepsilon,T) = U\bigl(\xi(\varepsilon)\bigr) + F_{\rm ph}\bigl(\xi(\varepsilon),T\bigr) + F_{\rm el}\bigl(\xi(\varepsilon),T\bigr),$$2 K for the majority of materials, and that near perfect agreement is obtained with five phonon spectra (Rostami et al., 2024). The detailed error summary states that E$$F_{\rm ZSISA}(\varepsilon,T) = U\bigl(\xi(\varepsilon)\bigr) + F_{\rm ph}\bigl(\xi(\varepsilon),T\bigr) + F_{\rm el}\bigl(\xi(\varepsilon),T\bigr),$$3Vib2 yields volumetric expansion and $$F_{\rm ZSISA}(\varepsilon,T) = U\bigl(\xi(\varepsilon)\bigr) + F_{\rm ph}\bigl(\xi(\varepsilon),T\bigr) + F_{\rm el}\bigl(\xi(\varepsilon),T\bigr),$$4 within $$F_{\rm ZSISA}(\varepsilon,T) = U\bigl(\xi(\varepsilon)\bigr) + F_{\rm ph}\bigl(\xi(\varepsilon),T\bigr) + F_{\rm el}\bigl(\xi(\varepsilon),T\bigr),$$5 of full QHA up to $$F_{\rm ZSISA}(\varepsilon,T) = U\bigl(\xi(\varepsilon)\bigr) + F_{\rm ph}\bigl(\xi(\varepsilon),T\bigr) + F_{\rm el}\bigl(\xi(\varepsilon),T\bigr),$$6 K for all but the softest or highly anharmonic cases, whereas E$$F_{\rm ZSISA}(\varepsilon,T) = U\bigl(\xi(\varepsilon)\bigr) + F_{\rm ph}\bigl(\xi(\varepsilon),T\bigr) + F_{\rm el}\bigl(\xi(\varepsilon),T\bigr),$$7Vib4 is essentially indistinguishable from QHA.

These developments clarify that ZSISA is not a single algorithm but a structural approximation that can be embedded in different numerical workflows. A plausible implication is that the main methodological frontier is no longer whether to use ZSISA at all, but how to parameterize $$F_{\rm ZSISA}(\varepsilon,T) = U\bigl(\xi(\varepsilon)\bigr) + F_{\rm ph}\bigl(\xi(\varepsilon),T\bigr) + F_{\rm el}\bigl(\xi(\varepsilon),T\bigr),$$8 efficiently enough that anisotropic QHA becomes routine for lower-symmetry systems.

5. Extension to polar solids and properties sensitive to internal coordinates

The 2026 ZnO study generalizes an approach originally developed for hcp metals to solids with both internal and external degrees of freedom, including non-vanishing piezoelectric and pyroelectric tensors (Gong et al., 5 Mar 2026). For wurtzite ZnO, the external parameters are $$F_{\rm ZSISA}(\varepsilon,T) = U\bigl(\xi(\varepsilon)\bigr) + F_{\rm ph}\bigl(\xi(\varepsilon),T\bigr) + F_{\rm el}\bigl(\xi(\varepsilon),T\bigr),$$9 and the internal parameter is the symmetry-allowed anion–cation shift $\xi(\varepsilon)$0 along the $\xi(\varepsilon)$1 axis. In this setting, ZSISA again defines

$\xi(\varepsilon)$2

where $\xi(\varepsilon)$3 is obtained by minimizing the static energy with respect to the internal coordinate at each external geometry. The Gibbs free energy under pressure is then minimized only with respect to the external parameters (Gong et al., 5 Mar 2026).

A central theoretical point is stated explicitly as a theorem attributed to Allan et al., J. Chem. Phys. 105, 8300 (1996): to linear order in $\xi(\varepsilon)$4, the external thermal expansion obtained by ZSISA is identical to the exact QHA result, but the internal thermal expansion is not (Gong et al., 5 Mar 2026). This sharply delineates the class of observables for which ZSISA is expected to perform well.

The ZnO study shows how this distinction propagates into derived tensors. The internal thermal expansion is defined as

$\xi(\varepsilon)$5

but under ZSISA the temperature dependence of $\xi(\varepsilon)$6 arises only indirectly through the external thermal strains because $\xi(\varepsilon)$7 does not respond directly to $\xi(\varepsilon)$8 at fixed $\xi(\varepsilon)$9 (Gong et al., 5 Mar 2026). Likewise, the piezoelectric tensor is decomposed into clamped-ion and internal-strain contributions,

$F_{\rm ph}$0

and the pyroelectric coefficient is written as

$F_{\rm ph}$1

In ZSISA, the first term in the pyroelectric coefficient vanishes because $F_{\rm ph}$2 has no direct $F_{\rm ph}$3-driven component; only the secondary term from lattice expansion survives (Gong et al., 5 Mar 2026).

Numerically, the study reports that ZSISA reproduces external thermal expansion $F_{\rm ph}$4 for ZnO to within a few percent, but underestimates the change of the internal parameter $F_{\rm ph}$5 by $F_{\rm ph}$6 at $F_{\rm ph}$7 K. As a consequence, the primary pyroelectric coefficient is significantly underestimated; the paper reports that $F_{\rm ph}$8 is $F_{\rm ph}$9 smaller in ZSISA and that the primary pyroelectricity is underestimated by about $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$00 because the direct vibronic shift of $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$01 is missed (Gong et al., 5 Mar 2026). By contrast, the piezoelectric coefficients $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$02 agree within $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$03 because $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$04 is only weakly temperature-dependent.

This establishes an important correction to a broad-brush view of ZSISA. It is not merely an accuracy-versus-cost trade-off in a generic sense; rather, it is systematically reliable for externally controlled thermoelastic observables and systematically less reliable for quantities controlled by the direct thermal evolution of internal coordinates.

6. Assumptions, failure modes, and methodological practice

The assumptions behind ZSISA are stated directly in the osmium study: ionic-position relaxations under phonon free-energy forces must be small, and the phonon contribution to the internal stress must be subdominant to the static DFT contribution (Gong et al., 21 Jul 2025). When these conditions hold, as reported for heavy transition-metal hcp systems, ZSISA reproduces FFEM to better than $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$05 while using typically an order of magnitude fewer phonon calculations; adding V-ZSISA can reduce the workload further by another factor of $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$06–$$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$07, at the expense of shape-coupling errors of order $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$08–$$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$09 (Gong et al., 21 Jul 2025).

The stated limitations are equally specific. The osmium and beryllium analyses both note that ZSISA may fail in very soft crystals, in materials with large internal relaxations, in systems with strong coupling to soft optical branches, or when genuine anharmonicity beyond QHA becomes important near high temperatures or structural phase transitions (Gong et al., 21 Jul 2025). The anisotropic ZSISA paper adds that, for strongly anharmonic systems and soft modes, beyond-QHA methods such as self-consistent phonons or molecular dynamics are needed above roughly the Debye temperature scale (Rostami et al., 12 Mar 2025).

The literature also gives explicit best-practice recommendations. For future QHA calculations on hcp metals, the osmium study proposes the following sequence: first compare frozen-ion and relaxed-ion elastic constants at $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$10 as a static-stress test; if static errors are $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$11–$$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$12, proceed with ZSISA; employ V-ZSISA by sampling phonons only along the $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$13 stress-pressure curve, with $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$14–$$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$15 points typically sufficient for hcp; validate the approximation with one FFEM calculation at the highest temperature of interest; use polynomial interpolation of degree $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$16–$$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$17 in volume and shape parameters; and document both ZSISA and V-ZSISA results alongside a single FFEM point (Gong et al., 21 Jul 2025).

A related practical lesson from the anisotropic study is that v-ZSISA should not be used when anisotropic expansion itself is the object of interest. That work reports, for example, that in wurtzite ZnO the two coefficients $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$18 and $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$19 become nearly identical in v-ZSISA-QHA, contrary to experiment, whereas the full anisotropic ZSISA–E$$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$20Vib2 restores $$F(\{C_\gamma\},T)=E_{\rm BO}(\{C_\gamma\})+F_{\rm vib}(\{C_\gamma\},T)$$21 (Rostami et al., 12 Mar 2025). This suggests a useful methodological distinction: v-ZSISA is a reduced model for volumetric thermodynamics, while full anisotropic ZSISA is the appropriate approximation for tensorial thermoelasticity.

Taken together, these studies position ZSISA as a controlled QHA approximation whose validity depends less on crystal symmetry alone than on the hierarchy of couplings among static restoring forces, vibrational stresses, and internal-coordinate relaxations. In systems where the phonon-induced internal forces are weak and the internal relaxation manifold is stiff, ZSISA is numerically close to FFEM. Where internal coordinates directly mediate the targeted observable, FFEM or an equivalent treatment of the internal free-energy landscape is required.