Harmonic Approximation: implementation and workflow

Purpose and scope

The harmonic-approximation workflow evaluates vibrational thermodynamic properties from phonon frequencies supplied at one or more fixed volumes. It implements the independent harmonic-oscillator expressions described in Harmonic Approximation, but it does not minimize the free energy with respect to volume. Every sampled volume remains an independent thermodynamic state.

This distinction is operationally important:

• HA answers what are the harmonic properties at each supplied volume?;

• QHA additionally asks which volume minimizes \(F(V,T)+PV\) at each pressure and temperature?

The HA workflow is therefore both a complete calculation in its own right and the first numerical stage of a QHA calculation.

Computational pipeline

The workflow follows the sequence

phonon YAML
    │
    ├─ validate volumes, static energies, frequencies, and q-point weights
    ├─ normalize q-point weights
    ├─ convert temperature to K and frequencies to Hz
    ├─ evaluate harmonic oscillator sums on the complete T × V grid
    ├─ convert properties to the requested output units
    └─ build ResultData → HDF5 / report / plots / table export

The numerical thermodynamic backend is vectorized with NumPy. Each property is evaluated over all temperatures, q-points, modes, and volumes without a Python loop over individual oscillators.

Scientific input contract

A normalized HA input contains at least:

volume

One value for every sampled structure, with shape (nvol,).

energy

Static electronic energy at every sampled volume, also with shape (nvol,).

frequencies

Mode-resolved phonon frequencies with shape (nq, 3 * natoms, nvol).

weights

One q-point weight for every q-point, with shape (nq,).

natoms and formula_units

Cell normalization information used by reports and molar conversions.

Quantas checks that the number of modes is exactly \(3N\), that every volume has one static energy and one complete phonon spectrum, and that the q-point and weight axes agree.

Multi-volume HA

HA is natively multi-volume. For every sampled volume \(V_i\), Quantas calculates

\[U_0(V_i),\quad U_{\mathrm{zp}}(V_i),\]

and the temperature-dependent arrays

\[U_{\mathrm{th}}(T,V_i),\quad U(T,V_i),\quad S(T,V_i),\quad F_{\mathrm{vib}}(T,V_i),\quad F(T,V_i),\quad C_V(T,V_i).\]

The result therefore preserves the full sampled temperature × volume surface. No volume is labelled as an equilibrium volume by HA.

q-point weights

The stored q-point weights are normalized internally by their sum before any oscillator sum is evaluated. Multiplying all weights by the same positive constant therefore does not change the result.

Quantas does not currently perform an independent symmetry analysis to infer q-point equivalence or multiplicity. Consequently:

• a full q-point set may use equal weights when scientifically appropriate;

• an irreducible set must carry the correct multiplicities from its source;

• weights imported from a code such as Phonopy are consumed as supplied;

• users should not invent multiplicities solely from visual similarity of q-point coordinates.

The format and provenance requirements are described in HA and QHA phonon YAML.

Treatment of zero and imaginary frequencies

The harmonic numerical functions exclude frequencies less than or equal to zero from the oscillator sums. This prevents undefined logarithms and Bose factors from contaminating all properties, but it must not be mistaken for a physical correction of an unstable structure.

A non-positive frequency may represent:

• a translational acoustic mode affected by numerical noise;

• an insufficiently converged force-constant calculation;

• an inconsistent supercell or sampling setup;

• a genuine dynamical instability.

Quantas does not decide which interpretation is correct. Before accepting a thermodynamic result, inspect how many modes were omitted, their magnitude, their q-point location, and whether the behavior persists across volumes. Removing a genuinely unstable mode can make a table finite without making the harmonic model scientifically valid.

Temperature grid and the zero-temperature limit

The public temperature range is defined by MIN MAX STEP. Quantas uses the same inclusive grid convention in HA and QHA:

np.arange(MIN, MAX + STEP, STEP)

Choose MAX consistently with STEP. When the interval is not an exact multiple of the step, this convention may include the next grid point beyond the nominal maximum.

At \(T=0\) the implementation uses explicit limits:

• \(U_{\mathrm{th}}=0\);

• \(S=0\);

• \(C_V=0\);

• \(F_{\mathrm{vib}}=U_{\mathrm{zp}}\).

The zero-point energy is evaluated once as a volume-dependent quantity and is stored with shape (1, nvol). Reports and plots broadcast it in temperature only for presentation; the native result is not duplicated.

Units and normalization

Before the oscillator sums are evaluated:

• temperatures are converted to kelvin;

• frequencies are converted to hertz;

• q-point weights are normalized.

The low-level harmonic functions return energy-like quantities in kJ mol^-1 and entropy/heat capacity in J mol^-1 K^-1. The workflow then converts them to the requested Quantas energy convention while preserving float64 values. Display precision is applied only by renderers.

Default behavior

HA defaults and rationale

Control	Default	Rationale
Temperature range	298.15 298.15 1.0 K	Produces one ambient-temperature state unless the user requests a thermodynamic curve.
Energy unit	Ha	Preserves the native electronic-structure energy scale used by the normalized input.
Volume length unit	A	Cell volumes are expressed as the cube of this length unit.
Frequency unit	cm^-1	Conventional output of vibrational calculations and spectroscopy.
Thermodynamic calculation	enabled	The public CLI always evaluates the requested HA properties.

The defaults define a safe minimal calculation, not a recommended production temperature range. The scientific question determines the appropriate lower and upper temperatures and resolution.

Performance and memory use

The dominant work and temporary-array size scale approximately as

\[N_T\,N_q\,N_{\mathrm{mode}}\,N_V.\]

The calculation evaluates several thermodynamic functions sequentially, so a large T × q × mode × volume product affects both execution time and peak memory.

Practical acceleration strategies are:

1. use a coarse temperature grid during input validation;

2. calculate only the temperature interval relevant to the study;

3. avoid creating many nearly redundant temperature points for a smooth plot;

4. verify the final grid only after the phonon dataset is scientifically sound;

5. use --benchmark to report backend timings when performance is being investigated.

Reducing the number of sampled volumes or q-points is not a neutral numerical optimization: it changes the physical input and must be justified at the phonon-calculation level.

Note

The current HA workflow has no numerical parameter with a default value of 512. Its resolution is controlled by the supplied phonon dataset and the requested temperature grid. A value named 512 encountered in a different Quantas module or in historical material must not be transferred to HA.

Reports, plots, and persistence

The calculation returns a generic quantas.api.common.ResultData envelope containing a typed HA payload. From the same result, the public API can build:

• neutral report tables;

• neutral plot specifications;

• native HDF5 output;

• selected property tables.

The CLI follows the same sequence and automatically writes a deterministic log for a significant calculation. Plotting occurs after the numerical run and does not alter or recompute the stored thermodynamic arrays.

Warnings and failure modes

The workflow stops before calculation when:

• required arrays are missing;

• volumes and static energies have inconsistent lengths;

• the phonon array does not have shape (nq, 3 * natoms, nvol);

• weights do not match the q-point axis;

• the sum of weights is not positive;

• the temperature range or step is invalid;

• a requested unit conversion is unsupported.

Non-positive modes are not a hard numerical failure because they are excluded from the sums. They remain a scientific validation responsibility, as explained above.

Decision guide

Practical HA decisions

Situation	Recommended action	Reason
First check of a new input	Run one or a few temperatures and inspect all volumes.	Detect shape, unit, and unstable-mode problems quickly.
Low-temperature study	Include \(T=0\) and use a finer grid where heat capacity changes rapidly.	The thermodynamic functions are most curved at low temperature.
High-temperature trend only	Use a wider but coarser preliminary grid, then refine only if needed.	Harmonic curves are usually smoother at high temperature.
Irreducible q-point input	Verify source-provided multiplicities before running.	Quantas normalizes weights but does not infer them.
Several non-positive optical modes	Stop and re-evaluate the phonon calculation.	Excluding them may conceal a real instability.
Need equilibrium volume or \(C_P\)	Continue with QHA.	HA does not minimize volume and directly provides only \(C_V\).

Related documentation

• Scientific theory: Harmonic Approximation

• Complete worked example: Harmonic Approximation

• Input specification: HA and QHA phonon YAML

• CLI syntax: quantas ha

• Public Python surface: Harmonic Approximation API