V-ZSISA: Volume-Constrained QHA for Thermoelasticity

• V-ZSISA is a volume-constrained form of the Zero Static Internal Stress Approximation that restricts calculations to the 0 K stress-pressure isotherm.
• It approximates thermoelastic properties by evaluating free energy and elastic constants along a one-dimensional path, significantly reducing computational costs.
• Empirical assessments in systems like hcp Be and osmium demonstrate that V-ZSISA delivers near-experimental accuracy while requiring far fewer strained configurations.

Volume-Constrained ZSISA (V-ZSISA, or v-ZSISA) is a volume-constrained form of the Zero Static Internal Stress Approximation used within the quasi-harmonic approximation (QHA) to reduce the cost of finite-temperature first-principles calculations. In the ab initio thermoelastic literature it is also described as the statically constrained quasi-harmonic approximation. Its defining feature is that equilibrium structures are not searched over the full multidimensional space of crystal parameters and internal coordinates; instead, the calculation is restricted to a reduced set of reference geometries tied to the equilibrium-volume path, typically the $0$ K hydrostatic stress-pressure isotherm, or, in the purely volumetric case, to a one-dimensional volume minimization in which the volume is optimized with the full free energy while the remaining degrees of freedom at fixed volume are optimized with the Born–Oppenheimer energy only (Gong et al., 2024, Rostami et al., 2024).

1. Definition and scope

In the QHA treatment of solids, the free energy depends on structural variables that may include lattice parameters, cell-shape ratios, and internal atomic coordinates. For hcp metals, the relevant external variables are commonly written as $\xi=(a,c/a)$; for more general crystals the variables may be partitioned into external parameters $\xi_i$ and internal parameters $u_k$. V-ZSISA is introduced precisely because a full QHA treatment over this multidimensional parameter space is expensive: each reference geometry may require phonon calculations and, in thermoelastic applications, multiple strained-configuration calculations as well (Gong et al., 2024, Gong et al., 21 Jul 2025).

The central approximation is geometric rather than purely variational. Standard ZSISA addresses how internal coordinates are treated under strain, whereas V-ZSISA further constrains how the equilibrium thermodynamic path is followed through structural-parameter space. In the hcp-metal implementation, one first determines the equilibrium crystal parameters at $0$ K under hydrostatic pressure, obtaining a one-dimensional set of structures along the “stress-pressure” isotherm, and then interpolates free energies or elastic constants only on that path rather than on the full two-dimensional $(a,c/a)$ grid. This reduced-reference-geometry strategy is motivated by the observation that, at least in the studied cases, the finite-$T$ isobar and the $0$ K stress-pressure path are often close enough that the resulting error is small (Gong et al., 2024, Gong et al., 21 Jul 2025).

2. Thermodynamic formalism

A general QHA formulation writes the Helmholtz free energy as

$$F(\{\lambda_A\},T)=U(\{\lambda_A\})+F_{vib}(\{\lambda_A\},T)+F_{el}(\{\lambda_A\},T),$$

where $\lambda_A$ collects the structural degrees of freedom, $\xi=(a,c/a)$0 is the static DFT energy, $\xi=(a,c/a)$1 is the vibrational free energy, and $\xi=(a,c/a)$2 is the electronic excitation term. In the generalized formalism for crystals with external and internal degrees of freedom,

$\xi=(a,c/a)$3

with $\xi=(a,c/a)$4. At finite pressure, the equilibrium structure is obtained by minimizing

$\xi=(a,c/a)$5

or, more generally under applied stress, the corresponding Gibbs functional with stress–strain coupling (Gong et al., 5 Mar 2026).

For volumetric thermal expansion, v-ZSISA specializes this general problem to a one-dimensional volume optimization. The pressure condition is written as

$\xi=(a,c/a)$6

so that the equilibrium volume satisfies

$\xi=(a,c/a)$7

This formulation makes explicit that the volume is optimized using the full free energy, while the other degrees of freedom at fixed volume are slaved to the Born–Oppenheimer minimum (Rostami et al., 2024).

In thermoelastic applications, the isothermal elastic constants are extracted from free-energy second derivatives,

$\xi=(a,c/a)$8

and, when the reference structure carries nonzero stress, they are converted to stress–strain elastic constants through the Barron–Klein correction. Along a hydrostatic stress-pressure path, this correction simplifies to

$\xi=(a,c/a)$9

Adiabatic constants are then obtained from

$\xi_i$0

These formulas are the basis for the V-ZSISA evaluation of finite-$\xi_i$1 elastic constants in the hcp-metal studies (Gong et al., 2024, Gong et al., 21 Jul 2025).

3. Relation to ZSISA and FFEM

V-ZSISA is not synonymous with ZSISA. ZSISA concerns the treatment of internal relaxations, whereas V-ZSISA adds a constraint on the equilibrium path or interpolation manifold used in finite-$\xi_i$2 calculations. FFEM, by contrast, removes the ZSISA simplification and minimizes the free energy directly with respect to internal coordinates at each strain and temperature (Gong et al., 2024, Gong et al., 21 Jul 2025, Gong et al., 5 Mar 2026).

Framework	Internal-coordinate treatment	Equilibrium-path treatment
ZSISA	Internal coordinates relaxed using static or total energy at each strain	No special volume-constrained path implied
V-ZSISA	ZSISA-like internal treatment plus reduced reference geometries	Properties evaluated along the stress-pressure curve or volume-consistent path
FFEM	Free energy minimized with respect to internal coordinates at each $\xi_i$3 and strain	No reduced-path approximation

For hcp Be, the distinction can be written explicitly for the strain $\xi_i$4. In ZSISA,

$\xi_i$5

so minimizing with respect to the internal displacement $\xi_i$6 gives

$\xi_i$7

and therefore

$\xi_i$8

FFEM replaces this by a free-energy minimization at each $\xi_i$9,

$u_k$0

yielding

$u_k$1

This is why FFEM is described as more general and more direct than either ZSISA or V-ZSISA (Gong et al., 2024).

The generalized ZnO study sharpens the distinction further. It states that, to lowest order—specifically, if $u_k$2 is a linear polynomial of $u_k$3—the external thermal expansion obtained from ZSISA is exact, but the internal thermal expansion is not correct. A common misconception is therefore that ZSISA-type approximations are uniformly accurate for all thermally induced structural changes; the cited results indicate that this is not so for internal-coordinate evolution, even when external lattice expansion is well reproduced (Gong et al., 5 Mar 2026).

4. Reduced-reference-geometry implementations

The practical form of V-ZSISA is a reduced-reference-geometry workflow. In the hcp Be implementation, the full interpolation used a $u_k$4 grid of crystal parameters $u_k$5, spanning pressures from about $u_k$6 kbar to $u_k$7 kbar, with phonons and $u_k$8 K elastic constants computed on all grid points. The practical V-ZSISA dataset, by contrast, consisted of only $u_k$9 selected points along the $0$0 K stress-pressure isotherm, with QHA elastic constants calculated in $0$1 of these geometries. Both the one-dimensional V-ZSISA interpolation and the full two-dimensional interpolation were fitted with a fourth-degree polynomial (Gong et al., 2024).

Component	Full-grid treatment	V-ZSISA treatment
Structural sampling	$0$2 grid in $0$3	11 points on the $0$4 K stress-pressure isotherm
QHA elastic-constant geometries	Grid-wide dataset	8 selected geometries
Interpolation	Fourth-degree polynomial on the 2D grid	Fourth-degree polynomial along the 1D path

The computational savings are substantial because each elastic-constant point is itself expensive. In the Be V-ZSISA-QHA workflow, each point required 5 strain types and 6 strain amplitudes each, giving $0$5 strained configurations, while FFEM required additional sampling of the internal coordinate $0$6 (Gong et al., 2024).

In the purely volumetric setting, v-ZSISA-QHA serves as a one-dimensional reference method for thermal expansion. The reference treatment typically used 7 phonon spectra spanning about $0$7 to $0$8 around $0$9, with more volumes in some cases such as +10% for Al, Cu, and Bi. A hierarchy of cheaper approximations was then introduced between the linear Grüneisen approach and v-ZSISA-QHA: E2Vib1 and E$(a,c/a)$0Vib1 require 2 phonon spectra, E$(a,c/a)$1Vib2 requires 3, and E$(a,c/a)$2Vib4 requires 5 (Rostami et al., 2024).

The same reduced-manifold idea appears in the generalized ZnO work, where figure captions explicitly describe thermoelastic properties as obtained within the V-ZSISA framework by interpolating results across five geometries along the stress-pressure curve. The paper does not re-derive a separate V-ZSISA theory, but its usage is explicitly that of a constrained or reduced-geometry ZSISA interpolation strategy rather than a full FFEM treatment (Gong et al., 5 Mar 2026).

5. Accuracy and empirical assessment

The available evidence presents V-ZSISA as a cost-reduction device whose error is often small but system-dependent. In hcp Be, QSA elastic constants obtained by interpolating only along the $(a,c/a)$3 K stress-pressure isotherm were found to be very close to those obtained from the full two-dimensional interpolation. For the internal-relaxation problem, ZSISA was also found to be a very good approximation to FFEM: at $(a,c/a)$4 K the reported deviations were $(a,c/a)$5 kbar ($(a,c/a)$6) and $(a,c/a)$7 kbar ($(a,c/a)$8). The same study emphasized that QHA is closer to experiment than QSA, because QSA includes thermal effects only through thermal expansion whereas QHA also includes the explicit temperature dependence of the free energy (Gong et al., 2024).

The osmium study is particularly informative because it makes V-ZSISA a stricter geometric test. In Be, the $(a,c/a)$9 kbar isobar is almost parallel to the $T$0 K stress-pressure curve; in Os, the $T$1 kbar isobar is almost perpendicular to the $T$2 K stress-pressure curve. Moreover, in osmium both $T$3 and $T$4 increase with temperature, whereas along the $T$5 K stress-pressure curve $T$6 increases while $T$7 decreases. Despite this markedly different geometry, the paper reports that the V-ZSISA effect is negligible on a large scale, though visible on enlarged plots. It also finds negligible deviations between ZSISA and FFEM, with maximum differences at $T$8 K of about 10 kbar, and concludes that QHA within V-ZSISA matches experiment better than QSA; indeed, QSA within V-ZSISA predicts the wrong sign for the temperature trend of $T$9 relative to experiment (Gong et al., 21 Jul 2025).

The volumetric thermal-expansion benchmarks on 12 materials provide a complementary assessment of the reference role of v-ZSISA-QHA. The linear Grüneisen approach is cheapest, requiring only two phonon spectra, but its accuracy deteriorates except at low temperatures. By contrast, the quadratic Taylor-based approximation E$0$0Vib2, which uses three phonon spectra, was found to determine thermal expansion in reasonable agreement with the v-ZSISA-QHA method up to $0$1 K for the majority of materials, while E$0$2Vib4, using five spectra, gave near-perfect agreement (Rostami et al., 2024).

6. Extensions, applicability, and limitations

The main modern extension of the V-ZSISA framework is to crystals with both external and internal degrees of freedom, including insulating polar solids. For wurtzite ZnO, the generalized QHA formalism uses external parameters $0$3 and $0$4 and one internal parameter $0$5. Within this setting, the paper compares ZSISA and FFEM for thermal expansion, elastic constants, dielectric constants, piezoelectric tensors, and pyroelectricity, while using V-ZSISA-style interpolation along selected stress-pressure geometries for thermoelastic properties (Gong et al., 5 Mar 2026).

This extension clarifies where V-ZSISA-type reductions are reliable and where they are not. The ZnO results state that external thermal expansion is fairly accurate in ZSISA, and many elastic constants and their qualitative temperature trends are reasonable, but the internal thermal expansion $0$6 is not correctly reproduced. FFEM is therefore better for internal-coordinate evolution and for derived pyroelectric properties. The finite-temperature pyroelectric coefficient is written as

$0$7

so any deficiency in the treatment of $0$8 propagates directly into the primary pyroelectric term (Gong et al., 5 Mar 2026).

The principal limitation of V-ZSISA in hcp metals is that it neglects the shift of $0$9 away from the $$F(\{\lambda_A\},T)=U(\{\lambda_A\})+F_{vib}(\{\lambda_A\},T)+F_{el}(\{\lambda_A\},T),$$0 K stress-pressure path as temperature changes. For Be this was judged acceptable, but the same work explicitly notes that this may not hold for other hcp metals. The osmium study, however, found that even when the geometry of the $$F(\{\lambda_A\},T)=U(\{\lambda_A\})+F_{vib}(\{\lambda_A\},T)+F_{el}(\{\lambda_A\},T),$$1 kbar isobar differs strongly from the $$F(\{\lambda_A\},T)=U(\{\lambda_A\})+F_{vib}(\{\lambda_A\},T)+F_{el}(\{\lambda_A\},T),$$2 K stress-pressure curve, the effect of V-ZSISA remains very small, and that a full-grid QHA interpolation requiring phonons at all distorted geometries is unnecessary at present for osmium (Gong et al., 2024, Gong et al., 21 Jul 2025).

A second limitation is computational scaling outside the reduced-path approximation. FFEM is feasible when the internal subspace is small, but the grid-point scaling becomes exponential as the number of internal degrees of freedom grows. V-ZSISA remains attractive precisely because it constrains the structural search to a low-dimensional path while preserving the main thermoelastic trends with much lower cost. In the current literature, its most defensible use is therefore as a controlled reduced-reference-geometry strategy: more accurate than crude volumetric surrogates, much cheaper than full multidimensional QHA or FFEM, and especially useful in high-pressure and high-temperature thermoelasticity where many state points must be sampled (Gong et al., 5 Mar 2026, Rostami et al., 2024).