Quasi-Harmonic Approximation: implementation and workflow

Purpose and scope

The QHA workflow converts a multi-volume phonon dataset into equilibrium and thermodynamic properties on a pressure-temperature grid. The physical model is described in Quasi-Harmonic Approximation; this page explains how Quantas represents the free-energy surface, how the available choices differ, and how to assess their numerical reliability.

QHA contains several legitimate routes because no single representation is optimal for every dataset. The main decisions are:

  • interpolate mode frequencies or integrated thermodynamic properties;

  • minimize a polynomial free-energy representation or an integrated equation of state;

  • derive \(\alpha_V\) from a mixed derivative, from mode Grüneisen parameters, or from numerical differentiation of \(V(P,T)\);

  • evaluate polynomial thermoelastic derivatives analytically or from a local volume grid.

A production calculation should treat these alternatives as sensitivity tests, not merely as interchangeable command-line settings.

Complete computational pipeline

multi-volume phonon YAML
    │
    ├─ validate V, U0, frequencies, weights, structure, and mode continuity
    ├─ inspect static E(V) and estimate the sampled pressure interval
    ├─ calculate HA properties on the sampled T × V grid
    ├─ fit one F(V,T) representation for every temperature
    ├─ evaluate the minimum of F(V,T) + P V for every requested pressure
    ├─ reconstruct thermodynamic properties at V(P,T)
    ├─ select alpha_V and calculate Cp, Ks, and macroscopic gamma
    ├─ reconstruct cell parameters and anisotropic expansion when available
    └─ persist arrays, diagnostics, warnings, and provenance in HDF5

The pressure-temperature result arrays use temperature as the first axis and pressure as the second axis.

Preflight inspection

Before a full run, use

quantas qha inspect input.yaml --eos BM3

The inspector fits the static energy-volume data using both a polynomial and, when requested, an integrated energy EOS. It reports the pressure associated with every sampled volume and the pressure range implied by each representation.

This preview is valuable because QHA reliability is controlled first by the volume support of the input. A dense pressure grid does not compensate for an equilibrium volume outside the sampled interval.

Inspect at least:

  • whether polynomial and EOS pressure estimates have the same trend;

  • whether either fit is classified as poor;

  • whether the intended pressure range is bracketed by the sampled volumes;

  • whether the reference volume lies inside a well-resolved energy basin;

  • whether the outermost volumes are sufficiently far from the expected minima.

Input requirements and mode continuity

QHA uses the same normalized phonon arrays as HA, but requires a genuine multi-volume series for interpolation and minimization. The input may also store a volume-constrained structural path.

The frequency scheme additionally depends on a mode-continuity status:

verified

Branch continuity across volume was explicitly checked.

assumed

The array order is accepted as continuous, but no explicit verification is recorded. Quantas permits the calculation and emits a warning.

unknown or unreliable

Frequency QHA is rejected because mode-by-mode volume fits would have no defensible correspondence.

The thermodynamic scheme does not require mode-by-mode continuity because it interpolates summed thermodynamic properties.

The sampled harmonic stage

For every requested temperature, Quantas first computes HA properties at all sampled volumes. The resulting total Helmholtz free energy

\[F(V_i,T)=U_0(V_i)+F_{\mathrm{vib}}(V_i,T)\]

is the common starting surface for both QHA schemes.

If harmonic thermodynamics cannot be evaluated, Quantas can continue with the static energy repeated at every temperature. This behavior preserves a useful static pressure-volume minimization and records a warning, but it is a degraded static-only calculation, not a complete QHA result. Temperature- dependent properties from such a run must not be interpreted as physical QHA predictions.

Frequency and thermodynamic schemes

freq: mode-resolved frequency interpolation

For every q-point and mode, Quantas fits

\[\nu_{qj}(V)\]

with a polynomial of degree frequency_degree. At each equilibrium volume, the frequencies are re-evaluated and all harmonic thermodynamic functions are recalculated from the oscillator expressions.

Advantages

  • preserves mode-resolved information;

  • enables mode Grüneisen parameters;

  • allows thermodynamic properties to be recalculated from one internally consistent interpolated spectrum;

  • for local polynomial derivatives, permits free energies to be regenerated around the equilibrium volume rather than read only from a global fit.

Limitations

  • requires reliable branch correspondence across volumes;

  • a mode crossing or accidental reordering can corrupt an otherwise smooth polynomial fit;

  • one fit is required for every q-point and mode;

  • fitted frequencies may become non-positive in extrapolation;

  • mode-resolved fitting and Grüneisen analysis increase time and memory use.

Use freq when the phonon branches are continuous and mode information is scientifically important.

td: thermodynamic-property interpolation

Quantas first calculates integrated HA properties at every sampled volume and then fits each property as a function of volume independently at every temperature. The fitted property is evaluated at the QHA equilibrium volume.

Advantages

  • does not require mode-by-mode branch continuity;

  • is insensitive to harmless permutations of phonon modes;

  • usually provides a robust route to integrated thermodynamic quantities;

  • avoids storing and differentiating mode-resolved fits.

Limitations

  • loses mode-resolved interpretation;

  • cannot provide mode Grüneisen parameters;

  • different thermodynamic properties are fitted separately, so consistency depends on the quality of all property fits;

  • high-order volume derivatives remain sensitive to interpolation choices.

Use td when branch tracking is unavailable or the scientific objective is limited to integrated thermodynamic properties.

Comparison

Frequency and thermodynamic QHA schemes

Aspect

freq

td

Interpolated object

Every mode frequency

Integrated HA properties

Mode continuity

Required

Not required

Mode Grüneisen parameters

Available

Unavailable

Main sensitivity

Branch ordering and frequency fits

Property-volume fits

Typical computational cost

Higher

Lower

Preferred use

Mode-resolved analysis

Robust integrated thermodynamics

Agreement between the two schemes for one material is evidence of robustness, not a universal guarantee. The method-comparison exercise in Quasi-Harmonic Approximation demonstrates the recommended check.

Representing the free-energy curve

At each temperature Quantas fits one free-energy model over the sampled volumes. The fitted model is reused for all requested pressures at that temperature. Increasing the number of pressure points therefore does not require refitting \(F(V,T)\), although later property evaluation still scales with the number of pressure states.

Polynomial minimization

The polynomial path fits \(F(V,T)\) in a centered and scaled volume coordinate. Scaling improves numerical conditioning without changing the physical volume derivatives.

At pressure \(P\), Quantas finds positive-curvature roots of

\[\frac{\partial F}{\partial V}+P=0.\]

If several local minima exist, the implementation selects the one connected to the sampled free-energy basin: it first identifies the sampled volume with the lowest \(F+PV\), then chooses the closest stationary minimum. Objective values break an exact distance tie.

Advantages

  • flexible local representation near the sampled basin;

  • no commitment to a specific EOS family;

  • efficient when the relevant minimum is well bracketed;

  • convenient for comparing alternative interpolation degrees.

Limitations

  • behavior outside the sampled interval is uncontrolled;

  • high polynomial degree can produce extra extrema or edge oscillations;

  • \(K'_T\) depends on a third volume derivative and is more sensitive than the equilibrium volume;

  • polynomial coefficients are numerical parameters, not material constants.

The default degree is three. The fit requires more volume points than fitted coefficients, and a scientifically useful fit normally needs additional redundancy beyond the mathematical minimum.

EOS minimization

The EOS path fits an integrated energy EOS to \(F(V,T)\) at every temperature and evaluates that fitted model at each pressure.

Advantages

  • imposes a physically structured compression curve;

  • gives interpretable parameters such as a reference volume and bulk modulus;

  • usually behaves more regularly over a wider compression interval;

  • supports covariance-based uncertainty propagation when the fit covariance is usable.

Limitations

  • results depend on EOS family and order;

  • parameters may be strongly correlated for a narrow volume interval;

  • a good residual norm does not guarantee reliable high-order derivatives;

  • EOS extrapolation remains extrapolation and must be flagged as such.

The default EOS tag is BM3. Alternative energy EOS forms are described in Equations of state and can be explored through the CLI reference.

Choosing polynomial or EOS minimization

Initial minimization choice

Dataset or objective

Start with

Why

Dense, symmetric sampling around the minimum

polynomial

Provides a flexible local representation.

Broad compression interval

EOS

A structured compression law is usually more stable.

Few sampled volumes

low-order model

Extra parameters are not independently resolved.

Strong disagreement between polynomial and EOS preview

neither blindly

Reinspect input range, convergence, and outliers.

Properties requested outside sampled volume support

extend the input

Changing the fit form cannot create missing physical information.

Polynomial degrees

The CLI exposes two principal degree controls.

--energy-degree

Sets the degree used for static energy, free energy, and thermodynamic- property volume fits. The default is three.

--frequency-degree

Sets the degree of every mode-resolved frequency-volume fit in the freq scheme. The default is three.

The public API also has a structural-path degree, default three, for the volume-dependent cell-shape reconstruction.

Increasing a degree may reduce residuals while worsening extrapolation and higher derivatives. Use the smallest degree that resolves reproducible curvature. Inspect residuals mode-by-mode or temperature-by-temperature, and compare the final physical properties rather than selecting a model from \(R^2\) alone.

Polynomial thermoelastic derivatives

After polynomial minimization, Quantas must calculate

\[K_T=V\frac{\partial^2F}{\partial V^2},\]

and

\[K'_T=-1-V \frac{\partial^3F/\partial V^3}{\partial^2F/\partial V^2}.\]

Two implementations are available.

local_grid — default

A small symmetric volume grid is built around every equilibrium volume. The default contains 5 points, with adjacent points separated by 0.05% of the central volume. The local free energies are fitted again to obtain the bulk properties.

With scheme=freq, local free energies are regenerated from the fitted phonon frequencies and static-energy model. With scheme=td, the local values come from the fitted global polynomial.

Advantages

  • evaluates curvature in the immediate neighborhood of the state;

  • provides a useful cross-check against derivatives of the global fit;

  • in the frequency scheme, uses the locally reconstructed oscillator free energy.

Limitations

  • adds local evaluations at every pressure-temperature state;

  • too small a separation can amplify floating-point or fit noise;

  • too large a separation no longer represents a local derivative;

  • \(K'_T\) remains intrinsically sensitive to the local polynomial degree.

analytic

The second and third derivatives are evaluated directly from the global free-energy polynomial.

Advantages

  • fastest route;

  • exact derivative of the selected global polynomial;

  • no additional local grid settings.

Limitations

  • inherits every high-order feature of the global polynomial;

  • can make \(K'_T\) especially sensitive to edge behavior or excessive degree.

Convergence test for local derivatives

At a small set of representative states, compare for example:

Suggested derivative convergence table

Points

Separation (%)

\(V\)

\(K_T\)

\(K'_T\)

3

0.05

…

…

…

5

0.05

…

…

…

7

0.05

…

…

…

5

0.025

…

…

…

5

0.10

…

…

…

analytic

—

…

…

…

The equilibrium volume should not change when only the derivative method is changed. Material changes in \(K_T\) or \(K'_T\) indicate that the curvature is not robustly resolved.

Note

There is no HA/QHA numerical control with a default value of 512 in the current source. The relevant QHA local derivative defaults are 5 points and 0.05% separation. The main global resolutions are the user-defined temperature and pressure grids. Do not tune an unrelated 512 value as though it controlled QHA convergence.

Three routes to volumetric thermal expansion

Quantas stores the available estimates separately and selects one requested method as the authoritative \(\alpha_V\). If the requested mixed or mode-Grüneisen value is unresolved at an individual state, the implementation uses the numerical volume derivative at that state when available. A source code array records which method supplied every value.

Mixed derivative — default

The implementation uses the Maxwell relation

\[\alpha_V =\frac{1}{K_T}\left(\frac{\partial S}{\partial V}\right)_T.\]

At every temperature, entropy on the sampled volume grid is fitted by a polynomial. Its analytical volume derivative is evaluated at every equilibrium volume and divided by \(K_T\).

Why it is the default

  • available for both freq and td;

  • based directly on the fitted free-energy thermodynamics;

  • avoids differentiating the final \(V(T)\) curve;

  • provides a smooth pressure-dependent estimate when the entropy-volume surface is well resolved.

Limitations

  • depends on the entropy-volume polynomial and its derivative;

  • inherits uncertainty in \(K_T\);

  • can become unreliable near volume boundaries or with too few volumes.

Mode-Grüneisen route

For scheme=freq, Quantas differentiates each fitted frequency polynomial:

\[\gamma_{qj}(V) =-\frac{V}{\nu_{qj}} \frac{\partial\nu_{qj}}{\partial V}.\]

The modes are weighted by their harmonic contribution to \(C_V\), and the thermal expansion is obtained from

\[\alpha_V=\frac{\bar\gamma C_V}{K_TV}.\]

Advantages

  • supplies mode-resolved interpretation;

  • identifies modes driving positive or negative expansion;

  • provides an independent thermodynamic consistency check.

Limitations

  • requires mode continuity and positive fitted frequencies;

  • non-positive or failed modes are excluded by default and recorded;

  • sensitive to derivatives of every frequency-volume fit;

  • more expensive than an integrated-property route.

Numerical volume derivative

The numerical route evaluates

\[\alpha_V=\frac{1}{V} \left(\frac{\partial V}{\partial T}\right)_P\]

from the final equilibrium-volume columns.

Advantages

  • intuitive and independent of entropy or mode derivatives;

  • useful as a cross-check;

  • remains available as the pointwise fallback for unresolved primary methods.

Limitations

  • requires at least two temperatures;

  • depends on temperature spacing and endpoint differentiation;

  • can amplify small irregularities in the minimized volume curve;

  • cannot supply a value for a single-temperature calculation.

Choosing an expansion method

Thermal-expansion decision guide

Situation

Primary choice

Recommended check

General QHA calculation

mixed derivative

Compare with numerical differentiation.

Verified continuous modes

mode Grüneisen

Compare weighted and macroscopic \(\gamma\).

No trustworthy mode ordering

mixed derivative

Use td and compare with numerical values.

Coarse temperature grid

mixed derivative

Numerical \(\alpha_V\) may be too step-dependent.

Disagreement among methods

no automatic winner

Inspect volume support, fit degrees, mode continuity, and derivatives.

Derived thermodynamic quantities

After selecting \(\alpha_V\), Quantas calculates

\[C_P-C_V=\alpha_V^2K_TVT,\]
\[C_P=C_V+(C_P-C_V),\]

and

\[K_S=K_T\frac{C_P}{C_V}.\]

At \(T=0\), and wherever \(C_V\) is too small for a stable ratio, Quantas sets \(K_S=K_T\).

The optional macroscopic Grüneisen parameter is

\[\gamma=\frac{\alpha_VK_TV}{C_V}.\]

The default low-heat-capacity threshold is 1% of the Dulong–Petit value \(3NR\). Below this threshold the ratio is reported as unresolved rather than allowing a numerically unstable value. At exactly zero temperature the stored macroscopic value is set to zero by convention.

Structural properties

When the input contains a volume-constrained structural series, Quantas builds a one-dimensional structural path in volume and evaluates it at \(V(P,T)\). Cell parameters and the deformation gradient are therefore reconstructed from static structures, while temperature enters through the QHA equilibrium volume.

Axial expansion is evaluated by the chain rule, for example

\[\alpha_a= \frac{\partial\ln a}{\partial\ln V}\,\alpha_V.\]

The tensor trace reproduces the selected volumetric expansion by construction. This is not a full anisotropic phonon QHA: it assumes that the cell-shape path is governed by the sampled static volume-constrained structures.

Volume support, extrapolation, and fit quality

Every equilibrium volume is classified as:

inside

More than 5% of the sampled interval width from either boundary.

near_boundary

Inside the interval but within 5% of one boundary.

outside

Outside the sampled volume interval.

The default API policy is to warn on extrapolation. A strict application can request failure instead. Near-boundary points deserve scrutiny because curvature and structural derivatives are already weakly constrained there.

Polynomial fits are classified as poor when, among other diagnostics, the design matrix is rank deficient, its condition number exceeds \(10^{12}\), or \(R^2<0.95\). EOS fits use the common Quantas fit quality contract. A converged fit with warnings is retained by the default warn quality policy; a stricter API workflow can stop on poor quality.

Failure policies and partial results

A local state can fail because:

  • the free-energy fit failed;

  • no physical minimum was found;

  • thermoelastic derivatives were invalid;

  • the selected EOS could not be evaluated at the requested pressure;

  • an equilibrium volume violated a strict extrapolation policy.

The CLI exposes three fit-failure policies:

continue

Record the failed state and proceed through the grid.

stop — default

Stop after the configured number of consecutive failures, default five. Previously completed states remain in the result.

raise

Raise immediately for applications that require all states to succeed.

A partially completed HDF5 result is useful diagnostic evidence, but failed or masked points must not be silently interpreted as a complete thermodynamic surface.

Uncertainty behavior

The default public options request covariance propagation where supported. In the EOS minimization path, the covariance of the temperature-specific EOS fit can be propagated to \(V\), \(K_T\), and \(K'_T\).

Advanced API users may request Monte Carlo propagation and control sample count, seed, confidence level, and the minimum accepted fraction of physical samples. Bootstrap propagation is not currently available for EOS pressure states; the workflow records a warning and omits those uncertainties rather than silently substituting another method.

Polynomial minimization does not provide an equivalent full parameter- covariance propagation for all final state quantities. Do not compare missing polynomial uncertainties with EOS uncertainties as though they represented the same statistical model.

Performance and practical acceleration

The total cost has several components.

Sampled harmonic thermodynamics

Scales approximately as \(N_TN_qN_{\mathrm{mode}}N_V\).

Frequency fitting

In freq, one volume polynomial is fitted for every q-point and mode.

Free-energy fitting

One polynomial or EOS is fitted for every temperature, not every pressure.

State evaluation

Scales with \(N_TN_P\); frequency thermodynamics and mode Grüneisen analysis add mode-resolved work at these states.

Structural reconstruction

Scales with the final pressure-temperature grid and is usually secondary.

A disciplined acceleration strategy is:

  1. run qha inspect before any large grid;

  2. test one or a few temperatures and pressures;

  3. disable mode Grüneisen analysis when it is not required;

  4. use td when mode-resolved information is unnecessary or unreliable;

  5. compare polynomial and EOS minimization on a reduced grid;

  6. validate local derivative settings at representative states;

  7. only then expand the final pressure-temperature domain.

Adding pressure points is often cheaper than adding temperature points because the fitted free-energy representation is shared across all pressures at one temperature. Nevertheless, every final state still requires property reconstruction, so very dense pressure grids are not free.

The following changes accelerate a run without changing the source phonon dataset:

  • coarser preliminary temperature and pressure steps;

  • --no-mode-gruneisen when mode analysis is not needed;

  • --no-gruneisen when the macroscopic ratio is not needed;

  • analytic polynomial derivatives for a controlled sensitivity test;

  • the td scheme for integrated properties.

Do not accelerate a production calculation by deleting volumes until the impact on fitted minima and derivatives has been quantified.

Programmatic validation

The public QHA API includes result-validation and comparison helpers. They can be used to compare two calculations on the same grid, for example freq against td or polynomial against BM3 minimization, without writing custom array-alignment code.

A method comparison should focus on physical quantities, not only fit metrics:

  • equilibrium volume;

  • \(K_T\) and \(K'_T\);

  • \(\alpha_V\) and its source code;

  • \(C_P-C_V\);

  • \(K_S\);

  • macroscopic and mode-weighted Grüneisen parameters;

  • structural extrapolation masks.