Harmonic Approximation

This tutorial calculates harmonic thermodynamic properties for periclase (MgO) over a series of eleven volumes. The same calculation is performed through the command-line interface and the public Python API.

The tutorial demonstrates an important property of the HA workflow: a single input may contain several volumes. Quantas evaluates the harmonic thermodynamic functions independently at every sampled volume; it does not select an equilibrium volume in this workflow. Equilibrium as a function of pressure and temperature is introduced later by the QHA tutorial.

Scientific objective

For each sampled volume and temperature, calculate:

• zero-point energy, \(U_{\mathrm{zp}}(V)\);

• thermal vibrational energy, \(U_{\mathrm{th}}(T,V)\);

• total internal energy, \(U(T,V)\);

• entropy, \(S(T,V)\);

• vibrational and total Helmholtz free energy;

• isochoric heat capacity, \(C_V(T,V)\).

Dataset

The example is derived from a CRYSTAL17 quasi-harmonic calculation for MgO. The normalized Quantas YAML contains:

• 2 atoms in the primitive cell;

• 1 formula unit per cell;

• 32 supplied q-points;

• 6 phonon branches at every q-point;

• 11 volumes from approximately 15.50 to 21.16 ų;

• one static energy and one phonon-frequency set for every volume.

The complete input is distributed at examples/qha/crystal-qha/mgo_b3lyp.yaml.

Download the MgO YAML input

A shortened excerpt is shown below. The complete file also contains the volume-dependent structures and all frequency arrays.

job: MgO Periclase (CRYSTAL17 QHA)
natom: 2
formula_units: 1
units:
  energy: Ha
  volume: angstrom^3
  frequency: cm^-1
  length: angstrom
mode_continuity: verified
mode_continuity_metadata:
  method: crystal-qha
  source: crystal
supercell:
  - [-2,  2,  2]
  - [ 2, -2,  2]
  - [ 2,  2, -2]
structure:
  representation: primitive
  reference_index: 6
  volume_series:
    volume: [15.495426716592775, 16.068019225910940, ..., 21.164485271468507]
qpoints: 32
phonon:
  - q-position: [0.0, 0.0, 0.0]
    weight: 1
    band:
      - frequency: [0.0, 0.0, 0.0, ...]
      # ... remaining branches and q-points ...

Note

Quantas uses the q-point weights supplied by the input. It does not derive symmetry multiplicities during an HA or QHA run. Input weights must therefore already represent the intended sampling convention.

The normalized tutorial input can be regenerated directly from the distributed native CRYSTAL QHA output:

quantas ha inpgen examples/qha/crystal-qha/mgo-b3lyp-crystal-qha.out \
    --interface crystal-qha \
    --jobname "MgO Periclase (CRYSTAL17 QHA)" \
    --output mgo_b3lyp.yaml

The generated file records CRYSTAL’s source-managed mode-continuity result. That metadata is not required by HA itself, but it allows the same normalized file to be reused by the QHA tutorial without weakening the input contract. See Phonon input generation and mode continuity for the complete generation procedure.

Temperature grid

The calculation uses temperatures from 0 to 1000 K in steps of 100 K. This coarse grid is chosen for a readable tutorial. Production calculations should use a resolution appropriate to the scientific question.

Running HA from the command line

From the repository root, run:

quantas ha run examples/qha/crystal-qha/mgo_b3lyp.yaml \
    --temperature 0 1000 100 \
    --output mgo_ha.hdf5 \
    --report mgo_ha.log \
    --verbosity standard \
    --no-progress \
    --force

The command produces:

mgo_ha.hdf5

Native scientific result containing normalized metadata and all HA arrays.

mgo_ha.log

Deterministic plain-text report.

The beginning of the report confirms that the full multi-volume dataset was accepted:

Input summary: 2 atoms, 32 q-points, 11 volumes

Selected options
----------------

           Option          |  Value
  -------------------------+--------
  Temperature minimum      | 0
  Temperature maximum      | 1000
  Temperature step         | 100
  Temperature unit         | K
  Energy unit              | Ha
  Volume unit              | A
  Frequency unit           | cm^-1
  Numerical precision      | float64

Understanding the HA report

Volume-only quantities

The static energy and zero-point energy depend on volume but not on the requested temperature grid. They are therefore reported once per volume:

Static and zero-point energies
------------------------------

     V      |          U0         |        Uzp
   (A^3)    |         (Ha)        |        (Ha)
------------+---------------------+-------------------
15.49542672 | -2.754467055837E+02 | 7.053334071905E-03
18.89686185 | -2.754628820004E+02 | 5.470199108461E-03
21.16448527 | -2.754585023425E+02 | 4.614423074024E-03

For this dataset, \(U_{\mathrm{zp}}\) decreases as the cell expands. The phonon spectrum softens with increasing volume, so the sum of \(h\nu/2\) terms becomes smaller.

Temperature-volume arrays

The report then identifies the stored shapes:

Key   Description                        Available  Shape
Uzp   Zero-point energy                  yes        (1, 11)
Uth   Thermal energy                     yes        (11, 11)
Utot  Internal energy                    yes        (11, 11)
S     Entropy                            yes        (11, 11)
Fvib  Vibrational Helmholtz free energy  yes        (11, 11)
F     Helmholtz free energy              yes        (11, 11)
Cv    Isochoric heat capacity            yes        (11, 11)

The first axis of the temperature-dependent arrays contains the 11 temperatures; the second contains the 11 sampled volumes. The zero-point energy retains its physically meaningful shape (1, 11) because it does not depend on temperature.

At the reference-like volume of 18.8969 ų, selected values are:

T (K)	\(U_{\mathrm{th}}\) (kJ mol⁻¹)	\(F_{\mathrm{vib}}\) (kJ mol⁻¹)	\(S\) (J mol⁻¹ K⁻¹)	\(C_V\) (J mol⁻¹ K⁻¹)
0	0.000	14.362	0.000	0.000
300	5.000	11.694	25.558	35.922
1000	36.238	-27.417	78.017	47.638

The zero-temperature vibrational free energy equals the zero-point energy. With increasing temperature, entropy lowers \(F_{\mathrm{vib}}\), while \(C_V\) approaches the classical high-temperature limit.

Exporting a property table

A selected property can be exported from the HDF5 result. For example:

quantas ha export mgo_ha.hdf5 \
    --property Cv \
    --unit J/mol \
    --output mgo_ha_Cv.dat

Download the complete HA heat-capacity table

The table preserves one column per sampled volume:

# Quantas HA table export
# Property: Cv - Isochoric heat capacity
# Units: J mol^-1 K^-1
#
T / K       V=15.495427 / A^3   ...   V=18.896862 / A^3   ...   V=21.164485 / A^3
0.00000000  0.000000000000E+00 ...   0.000000000000E+00 ...   0.000000000000E+00
300.000000  2.980199098224E+01 ...   3.592223697777E+01 ...   3.903139950331E+01
1000.00000  4.670869809859E+01 ...   4.763783551464E+01 ...   4.802494090781E+01

Generating selected HA plots

Isochoric heat capacity

quantas ha plot mgo_ha.hdf5 \
    --property Cv \
    --unit J/mol \
    --output figures/mgo_ha_cv \
    --preset publication \
    --format png \
    --dpi 180

Harmonic \(C_V(T,V)\) for all sampled volumes. Larger volumes have softer frequencies and therefore reach a given fraction of the high-temperature limit at lower temperature. At high temperature the curves converge toward the same classical limit.

Vibrational Helmholtz free energy

quantas ha plot mgo_ha.hdf5 \
    --property Fvib \
    --unit kJ/mol \
    --output figures/mgo_ha_fvib \
    --preset publication \
    --format png \
    --dpi 180

Vibrational Helmholtz free energy as a function of temperature. Expanded cells have lower-frequency phonons and acquire a more negative entropic contribution. This changing volume dependence is the physical ingredient used by the quasi-harmonic approximation.

Volume sections and V–T maps

When at least two sampled volumes are available, the same stored HA arrays can be viewed along the complementary natural variable without recomputing the harmonic sums. The selected temperatures must be exact points of the stored grid:

quantas ha plot mgo_ha.hdf5 \
    --property F \
    --axis volume \
    --temperature 0 \
    --temperature 300 \
    --temperature 1000 \
    --2d \
    --output figures/mgo_ha_free_energy_sections

The line plot reports \(F(V)\) at the selected temperatures. The optional contour plot reports \(F(T,V)\) on the complete native grid. Quantas does not interpolate missing volumes or temperatures during plot construction.

Running the same calculation from Python

The complete API script is distributed with the examples:

Download the HA API script

"""Reproduce the HA MgO tutorial through the public Quantas API."""

from __future__ import annotations

import argparse
from pathlib import Path

from quantas.api import ha, rendering


def main() -> None:
    """Run HA, write HDF5/report files, and render selected figures."""
    parser = argparse.ArgumentParser()
    parser.add_argument("input", type=Path)
    parser.add_argument("--output-dir", type=Path, default=Path("ha_tutorial"))
    args = parser.parse_args()

    args.output_dir.mkdir(parents=True, exist_ok=True)

    options = ha.Options(
        temperature_min=0.0,
        temperature_max=1000.0,
        temperature_step=100.0,
        energy_unit="Ha",
        volume_unit="A",
        frequency_unit="cm^-1",
        temperature_unit="K",
    )

    result = ha.run(args.input, options=options)
    payload = ha.get_result(result)

    report = rendering.render_tables(ha.build_report(result))
    report_path = args.output_dir / "mgo_ha_api.log"
    report_path.write_text(report, encoding="utf-8")

    hdf5_path = ha.write_result(
        result,
        args.output_dir / "mgo_ha_api.hdf5",
        report_text=report,
    )

    plot_requests = (
        ("Cv", "J/mol", "mgo_ha_cv_"),
        ("Fvib", "kJ/mol", "mgo_ha_fvib_"),
    )
    for property_name, unit, prefix in plot_requests:
        collection = ha.build_plots(
            result,
            properties=property_name,
            unit=unit,
        )
        rendering.render_plots(
            collection,
            output_dir=args.output_dir,
            filename_prefix=prefix,
            preset="publication",
            dpi=180,
            close=True,
        )

    print(f"Temperatures: {payload.temperature.size}")
    print(f"Volumes: {payload.volume.size}")
    print(f"Cv shape: {payload.isochoric_heat_capacity.shape}")
    print(f"Written: {hdf5_path}")
    print(f"Written: {report_path}")


if __name__ == "__main__":
    main()

Run it with:

python examples/ha/tutorial_api.py \
    examples/qha/crystal-qha/mgo_b3lyp.yaml \
    --output-dir ha_tutorial

The relevant public steps are:

1. construct quantas.api.ha.Options;

2. call quantas.api.ha.run();

3. obtain the typed payload with quantas.api.ha.get_result();

4. build deterministic tables with quantas.api.ha.build_report();

5. write HDF5 with quantas.api.ha.write_result();

6. build neutral plots with quantas.api.ha.build_plots();

7. render them through quantas.api.rendering.

The script reports:

Temperatures: 11
Volumes: 11
Cv shape: (11, 11)
Written: ha_tutorial/mgo_ha_api.hdf5
Written: ha_tutorial/mgo_ha_api.log

Scientific interpretation

This HA calculation provides thermodynamic functions at fixed volumes. It shows that the volume dependence of the vibrational spectrum is substantial, but HA alone does not determine which volume is stable at a given pressure and temperature. The next tutorial combines the static and vibrational free energies, minimizes \(F(V,T)+PV\), and obtains equilibrium QHA properties.

Reproducibility checkpoints

For the distributed input and the options used above:

• the input contains 11 volumes and 32 supplied q-points;

• Uzp has shape (1, 11);

• every temperature-dependent property has shape (11, 11);

• at 18.89686185 ų, \(U_{\mathrm{zp}}\) is approximately 14.3620 kJ mol⁻¹;

• at 300 K and the same volume, \(C_V\) is approximately 35.9222 J mol⁻¹ K⁻¹.