Anisotropic thermal expansion from free-energy minimization#

This page computes the anisotropic thermal expansion of a crystal, with one expansion coefficient per crystal axis rather than one for the volume. The ordinary quasi-harmonic approximation (QHA) minimizes the free energy along a single path in volume. This page samples the lattice parameters on a grid and minimizes the free energy over them directly, so the axes are free to expand by different amounts.

Steps 0 to 4 are a QHA calculation. The phonons come from displaced supercells computed with the calculator, and no machine-learning potential (MLP) is involved. Most calculations need only these steps.

Step 5 replaces the harmonic free energies, for a crystal whose anharmonicity the harmonic approximation cannot carry. An MLP is trained at each grid point and gives force constants that change with temperature. The free energies from those force constants enter the same minimization in place of the force sets, so this step is no longer a QHA calculation.

The details of the implementation are described in https://arxiv.org/abs/2609.24336.

$a$, $b$ and $c$ on this page are the lattice parameters of the standardized conventional unit cell, never of the primitive cell. Step 0 explains why the recipe is built on that cell.

The free lattice degrees of freedom (DOF) are detected from the symmetry: one for cubic ($a$), two for hexagonal, tetragonal and rhombohedral ($a, c$), and three for orthorhombic ($a, b, c$). Cell angles are held fixed, so monoclinic and triclinic crystals are out of scope. This page uses $(a, c)$ throughout as a concrete example; substitute the free DOF of your system. The lattice parameters and axial thermal expansions are produced for any of the supported systems. Contour maps of the free energy $F(a, c)$ at a fixed temperature – the surface whose minimum gives those lattice parameters – are drawn only when there are exactly two free DOF.

Warning

This workflow is experimental. Everything on this page works and is tested, but the interfaces are not settled. The command-line options, the aniso_qha_dataset.hdf5 layout, and the phonopy.qha.anisotropic and phonopy.qha.anisotropic_dataset APIs may change in a backward-incompatible way between releases, without a deprecation period; options have already been added and removed as the recipe was used on real systems.

So rebuild the dataset from the calculator outputs rather than relying on an old file being readable, keep beside every result the commands that produced it, and pin the phonopy version whenever a campaign has to stay reproducible across releases.

Four commands do the work, one per step:

phonopy-init prepares the equilibrium reference of step 0, and collects the forces of the phonon grid in step 2.
phonopy-strain-cells samples the strained cells of step 1.
phonopy-anisotropic-qha-dataset gathers the calculator outputs into the intermediate dataset of step 3.
phonopy-anisotropic-qha runs the analysis of step 4.

Step 4 runs the analysis in one command. The API script printed under it does the same thing with more control.

Prerequisites: h5py, symfc, a VASP setup, and pypolymlp for the variant of step 5.

All lengths are in the native length unit of the input cell (Angstrom for VASP); no unit conversion is applied by the tools.

Note

VASP is the only calculator this workflow runs with. The commands and the helper scripts assume VASP inputs and outputs (POSCAR, vasprun.xml, vaspout.h5).

The binding constraint is the static grid. The internal energy $U(a_i, c_i)$ and the electronic states are read from VASP outputs, and phonopy has no interface yet for reading the static single-point energy of the other calculators. So phonopy-anisotropic-qha-dataset stops early on a reference naming one of them, rather than building a dataset with $U(a_i, c_i)$ missing. Nothing else stands in the way: the readers for the other calculators are simply not written yet.

The free energy#

The free energy minimized at each temperature is

\[F(a, c; T) = U(a, c) + F_\mathrm{ph}(a, c; T) + F_\mathrm{el}(a, c; T),\]

where the electronic term $F_\mathrm{el}$ is optional. $U(a, c)$ is the static internal energy, computed by the calculator on the static grid of step 1. $F_\mathrm{ph}(a, c; T)$ is the phonon free energy of the harmonic crystal,

\[F_\mathrm{ph}(a, c; T) = \sum_{\mathbf{q}\nu} \left[ \frac{\hbar \omega_{\mathbf{q}\nu}}{2} + k_\mathrm{B} T \ln \left( 1 - e^{-\hbar \omega_{\mathbf{q}\nu} / k_\mathrm{B} T} \right) \right],\]

where the frequencies $\omega_{\mathbf{q}\nu}$ come from that grid point’s force constants of step 2, and the sum runs over the modes at the q points of the sampling mesh. Steps 0 to 4 use this expression as written. Step 5 replaces it with the SSCHA free energy.

$a$ and $c$ are continuous above: $F$ is written for any lattice. The calculator returns values only at the sampled cells $(a_i, c_i)$ of step 1. Step 4 therefore fits a surface through those values and minimizes the surface, and that is what makes $a(T)$ and $c(T)$ continuous functions of temperature.

Overview#

The boxes are jobs run by phonopy tools or the API, the hexagons are calculator runs, and the rounded nodes are input and intermediate data. This is the recipe of steps 0 to 4, with the phonons from the calculator:

        flowchart TD
    EQ(["Equilibrium cell<br/>(phonopy_disp.yaml)"])
    EQ --> SC["phonopy-strain-cells<br/>(a, c grid)"]
    SC --> RELAX{{"calculator relax + static"}}
    RELAX --> SGRID(["static-grid/grid-NNN<br/>U, F_el"])

    SGRID --> PD["generate displacements<br/>per relaxed cell"]
    PD --> CALCF{{"calculator forces"}}
    CALCF --> PGRID(["phonon-grid/grid-NNN<br/>disp-*"])

    SGRID --> BUILD["phonopy-anisotropic-qha-dataset<br/>--static ... --phonon ..."]
    PGRID --> BUILD
    BUILD --> DS(["aniso_qha_dataset.hdf5<br/>cells, U, F_el,<br/>displacements, forces"])

    DS --> ANA["phonopy-anisotropic-qha"]
    ANA --> RES(["a(T), c(T),<br/>alpha_a, alpha_c,<br/>F(a,c) maps"])

0. The equilibrium reference (`phonopy_disp.yaml`)#

Step 0 is done once, before the grid exists; steps 1 to 4 are the calculation itself. Start from one relaxed equilibrium cell and turn it into the reference phonopy_disp.yaml with phonopy-init. Every later step reads this file:

% phonopy-init -c REFERENCE_UNITCELL -d --dim 4 4 4

REFERENCE_UNITCELL must be the standardized conventional cell, whose lattice vectors are the crystal axes a, b and c in that row order.

The conventional cell is what makes the grid small. The number of grid points is set by the number of free lattice DOF, and the symmetry of the conventional cell is what reduces that number: two for hexagonal instead of three, one for cubic instead of three. The primitive cell of a centred lattice hides that symmetry in its rows, and the same crystal would then be sampled over more DOF than it has.

The free lattice DOF are taken per row, so a primitive cell of a centred lattice cannot be used at all: its rows are centring vectors rather than crystal axes. For body-centred tetragonal, for example, all three primitive rows have the same length, and scaling them would only change the volume, never $c/a$. A rhombohedral cell must likewise be given in the hexagonal setting. phonopy-strain-cells rejects such a cell; if in doubt, take the BPOSCAR written by phonopy-init –symmetry. That file is the conventional cell. A conventional cell that is merely rotated in Cartesian space is fine.

Everything the recipe writes follows from this choice. phonopy-strain-cells strains the unit cell of phonopy_disp.yaml and writes unit cells; the displaced supercells are built on them; and the lattice parameters reported at the end are theirs. The primitive cell enters only as the normalization of the energies, which “Run the anisotropic QHA” covers.

--dim fixes the supercell matrix, which phonopy-init records together with the unit cell, the primitive matrix and the calculator. phonopy-strain-cells reads the equilibrium cell and calculator from it; phonopy-anisotropic-qha-dataset reads the calculator, and takes the free lattice DOF from the symmetry of that cell (which lengths are independent – $a, c$ with $b = a$ for hexagonal). Keep --dim consistent with the phonon-grid supercell in step 2.

1. Build the static grid (internal energy U)#

Sample strained unit cells over the free lattice DOF, then relax and run a static single point for each with the calculator. The static grid supplies $U(a_i, c_i)$ and, optionally, the electronic states for $F_\mathrm{el}$.

# Inspect the free lattice DOF first (no ranges -> DOF report):
% phonopy-strain-cells phonopy_disp.yaml

# --grid N is the number of grid points per free axis (5 -> 5 x 5 = 25 cells);
# one N per free DOF gives a rectangular grid, e.g. --grid 5 6 -> 30 cells.
% phonopy-strain-cells phonopy_disp.yaml --a 3.168 3.232 --c 5.148 5.252 --grid 5
Wrote 25 strained unit cell(s) as unitcell-001 .. unitcell-025 in vasp format.
Grid sampling: 5 x 5 over (a, c).
  Main diagonal (5 cells), the --compare-eos volume path:
    a  c   c/a
    3.1680  5.1480   1.6250
    3.1840  5.1740   1.6250
    3.2000  5.2000   1.6250
    3.2160  5.2260   1.6250
    3.2320  5.2520   1.6250
Provenance written to strain_cells.yaml

The main diagonal of the grid, printed above with each cell’s $c/a$, is the volume path that phonopy-anisotropic-qha --compare-eos fits. Equal fractional ranges and equal counts keep $c/a$ constant along it, which is the cleanest input to the cross-check. Unequal ranges or counts still give a path, but $c/a$ then varies along it, and the fit mixes a change of shape into what it reads as a change of volume.

The cells are written as strained unit cells unitcell-NNN. When primitive_matrix is not the identity, the primitive cell is a different cell, and the same run writes it as primcell-NNN too:

Wrote 25 strained unit cell(s) as unitcell-001 .. unitcell-025 in vasp format.
Also wrote the primitive cell of each as primcell-001 .. primcell-025
(primitive 2 atoms, unit cell 8 atoms). The static single point may be run on
either cell.

The test is on the matrix, not on the atom counts, so a primitive cell taken in an unusual setting is written as well even though it holds the same atoms. The hexagonal example above is the other case: auto resolves to the identity there, the two cells are one cell, and only the unit cells are written.

The two files are one grid point. Run the static single point on whichever suits the calculation and pass that one to the builder; “Run the anisotropic QHA” says how $U$ is put on one normalization afterwards.

The primitive cells take their atoms from phonopy_disp.yaml rather than from a fresh symmetry search. The strain scales the crystal axes and leaves the fractional coordinates alone, so only the lattice is rebuilt. The atom order is then the same at every grid point, and does not move with the symmetry tolerance.

Each run also writes strain_cells.yaml, recording the ranges, the grid shape and the free-DOF lengths of every cell, with the primcell-NNN name beside each unitcell-NNN when both were written. Nothing in the run is random, so the same command reproduces the same cells.

For each grid point, on whichever of the two cells you chose:

Relax the internal coordinates if the structure has free internal parameters (e.g. the wurtzite u). A crystal with no internal DOF (e.g. HCP) skips this step – the strained cell is already the relaxed cell.
Run a static single point. Three settings matter.
- ISIF >= 2, if you also want the stress.
- Write vaspout.h5, if you want the electronic states for $F_\mathrm{el}$. A run that writes only vasprun.xml carries no eigenvalues, and the dataset is then built without $F_\mathrm{el}$.
- Sample the Brillouin zone on a Gamma-centred regular mesh. The k points and the mesh are then stored with the states, and $F_\mathrm{el}$ is integrated by the linear tetrahedron method rather than summed over k points.
The mesh a static single point already needs is dense enough for the tetrahedron method, so no extra one has to be chosen. Where a denser mesh is wanted for the electronic states alone, add a KPOINTS_OPT block: the states are read from that block in preference to the SCF mesh.
Place the output in static-grid/grid-NNN/ (one directory per grid point, containing vaspout.h5 or vasprun.xml). Any layout works – the builder is given the paths explicitly – but a name that sorts in the sampling order keeps the two grids easy to pass together, and no index file is needed. Sorting in that order also matters for --compare-eos: the builder recognizes a tensor grid only when the cells reach it in the order they were generated, and records its shape for the main-diagonal volume path.

Script 1 below lays out those directories. Edit the paths at the top of it and run it. It writes the POSCARs only; the rest of the VASP input – INCAR, POTCAR, KPOINTS – has to be distributed separately.

Script 1 – the static-grid input POSCARs, from the strained cells above#

import glob
from pathlib import Path
from phonopy.interface.calculator import read_crystal_structure, write_crystal_structure

# Strained cells from phonopy-strain-cells. Put "primcell-*" here instead to
# run the static grid on the primitive cells.
CELLS = "unitcell-*"
STATIC_GRID = "static-grid"

for path in sorted(glob.glob(CELLS)):
    idx = int(Path(path).stem.split("-")[-1])
    cell, _ = read_crystal_structure(path, interface_mode="vasp")
    static_dir = Path(STATIC_GRID) / f"grid-{idx:03d}"
    static_dir.mkdir(parents=True, exist_ok=True)
    write_crystal_structure(static_dir / "POSCAR", cell, interface_mode="vasp")
    print(f"grid-{idx:03d}: static POSCAR")

Then relax (if the crystal has internal DOF) and run the static single point in each static-grid/grid-NNN/.

2. Compute the phonons (phonon grid)#

For each relaxed static-grid cell, generate displaced supercells and compute their forces with the calculator.

Place the results as phonon-grid/grid-NNN/ each containing the phonopy_disp.yaml and the per-displacement subdirectories disp-001/, disp-002/, … (each with vaspout.h5 or vasprun.xml). The names need not match those under static-grid/, since the two grids are paired by the order they are passed in, but matching names make that order easy to get right.

Script 2 below reads each static-grid/grid-NNN/CONTCAR, the relaxed structure, which equals the input POSCAR when there is no internal DOF. So run it only after the static grid is done. Edit the paths at the top of it, and distribute the rest of the VASP input separately, as in step 1.

SUPERCELL_MATRIX and DISTANCE must both be the same at every grid point, for one reason. The analysis fits a surface through the $F_i(T)$ of the grid and differentiates it. Either setting biases $F_\mathrm{ph}$ a little – the supercell size through how far the force constants reach, the displacement distance through the anharmonic error of the finite difference. A bias that is equal everywhere shifts the surface without tilting it, and the axial expansions do not see it. A bias that changes from one grid point to the next tilts the surface, and the axial expansions follow the tilt.

So keep SUPERCELL_MATRIX at the --dim of step 0, which is also the value stored when a dataset is built from the static grid alone. DISTANCE is the displacement distance in Angstrom, 0.03 here against phonopy’s own default of 0.01. Nothing checks either of them: the builder reads the supercell matrix from each phonon grid point separately, and never compares the grid points with one another or with step 0.

Script 2 – the phonon grid, from the relaxed static-grid cells#

from pathlib import Path
import phonopy
from phonopy.interface.calculator import read_crystal_structure, write_crystal_structure

STATIC_GRID = "static-grid"  # relaxed cells at static-grid/grid-NNN/CONTCAR
PHONON_GRID = "phonon-grid"
N_GRID = 25  # cells written by phonopy-strain-cells in step 1
SUPERCELL_MATRIX = [4, 4, 4]
DISTANCE = 0.03

for idx in range(1, N_GRID + 1):
    contcar = Path(STATIC_GRID) / f"grid-{idx:03d}" / "CONTCAR"
    cell, _ = read_crystal_structure(str(contcar), interface_mode="vasp")
    ph = phonopy.Phonopy(cell, supercell_matrix=SUPERCELL_MATRIX, calculator="vasp")
    ph.generate_displacements(distance=DISTANCE)
    phonon_dir = Path(PHONON_GRID) / f"grid-{idx:03d}"
    phonon_dir.mkdir(parents=True, exist_ok=True)
    ph.save(phonon_dir / "phonopy_disp.yaml")

    for k, sc in enumerate(ph.supercells_with_displacements, 1):
        disp_dir = phonon_dir / f"disp-{k:03d}"
        disp_dir.mkdir(parents=True, exist_ok=True)
        write_crystal_structure(disp_dir / "POSCAR", sc, interface_mode="vasp")
    print(f"grid-{idx:03d}: {len(ph.supercells_with_displacements)} disp")

Then build the intermediate dataset, which is step 3.

3. Build the intermediate dataset#

The analysis of step 4 reads aniso_qha_dataset.hdf5, and so does step 5. The grid points are given as two path lists, --static and --phonon, and paired by position after shell expansion:

Note

The builder reads VASP outputs only, and stops early on a reference naming another calculator. The opening note of this page says why.

% phonopy-anisotropic-qha-dataset phonopy_disp.yaml \
    --static static-grid/grid-{001..025}/ \
    --phonon phonon-grid/grid-{001..025}/ \
    -o aniso_qha_dataset.hdf5
Reading pre-computed forces for 25 grid point(s)
  grid 1 U=... eV n_disp=...
  ...
  grid 25 U=... eV n_disp=...
Grid shape [5, 5] recorded for the main-diagonal path.
Wrote 25 grid point(s) to aniso_qha_dataset.hdf5

The entries may be files rather than directories. Use that form when the calculations were laid out by something other than Scripts 1 and 2:

% phonopy-anisotropic-qha-dataset phonopy_disp.yaml \
    --static runs/*/static/vaspout.h5 \
    --phonon runs/*/phonons/phonopy_params.yaml \
    -o aniso_qha_dataset.hdf5

Any names work, and one list can mix directories and files. The --static and --phonon lists must have the same number of entries.

The grid points must have the same primitive matrix, supercell matrix and chemical symbols.

build_aniso_qha_dataset builds the same dataset from objects rather than from directories. Give it one Phonopy per grid point, each carrying that point’s cell and its displacements and forces, and the static internal energy of each point:

Building the dataset through the API#

from phonopy.qha.anisotropic_dataset import (
    build_aniso_qha_dataset,
    write_aniso_qha_dataset,
)

dataset = build_aniso_qha_dataset(
    phonopys,
    internal_energies,  # eV per primitive cell
    electronic_structures=electronic_structures,  # optional, for F_el
)
write_aniso_qha_dataset(dataset, "aniso_qha_dataset.hdf5")

Pass the internal energies per primitive cell, so that they match the phonon calculation.

What the builder reads#

For each grid point the builder reads:

the static single point, giving the internal energy $U(a_i, c_i)$. It may have been run on the grid point’s unit cell or on its primitive cell, and the builder accepts either; anything else is reported as a mis-pairing. With --phonon given, the grid-point cell is the phonon entry’s cell rather than this one, so a relaxation carried into the phonon grid is honored. Without --phonon this cell becomes the grid point itself, and it then has to be the conventional unit cell, since the free lattice DOF are read from its rows. The electronic states for $F_\mathrm{el}$ are read automatically from the same vaspout.h5 when it carries the eigenvalues (a static point written with only vasprun.xml is built without $F_\mathrm{el}$; pass --no-electronic to skip them deliberately). A directory entry is resolved to the VASP output it holds, and vaspout.h5 is used in preference to vasprun.xml.
the phonon grid point, in one of two forms. A directory holding phonopy_disp.yaml and the per-displacement disp-* subdirectories: the builder reads each disp-* calculator output itself, so no FORCE_SETS or phonopy_params.yaml is needed. Or a phonopy.yaml-like file carrying forces that phonopy-init -f has already collected. This second form is the simpler route when the calculations were not laid out by Script 2:
```
% phonopy-init --sp -f disp-*/vasprun.xml   # -> phonopy_params.yaml
% phonopy-init -f disp-*/vasprun.xml        # -> FORCE_SETS, beside phonopy_disp.yaml
```
Pass the resulting phonopy_params.yaml. A phonopy_disp.yaml with its FORCE_SETS beside it works too (–sp merges the two into one file). A file with no forces and no neighboring FORCE_SETS is rejected rather than silently producing an empty grid point. Either form supplies the per-grid-point supercell / primitive matrices.

The positional phonopy_disp.yaml is the equilibrium reference; it supplies the free lattice DOF metadata and the calculator. The grid-point index recorded in the dataset is the position in the list, a label only, since the analysis reads the lattice parameters from each stored cell.

How the ordering is checked#

The builder pairs --static and --phonon by position after shell expansion, and takes the disp-* subdirectories of one grid point in sorted order. It checks both against the structures in the calculator outputs.

Across grid points, the lattice of each static single point must match the cell of the phonon grid point it is paired with. A mismatch stops the command and names both paths. Two mistakes are caught this way.

A grid point missing from each list. The list lengths still match, so nothing else notices, and the $U$ of one lattice would be combined with the forces of another.
A static single point run on a supercell. Its $U$ would then be on the wrong normalization.

Within one grid point, sorted order is only a guess at the displacement order: disp-1, disp-10, disp-2 sorts differently from how it counts. Each calculator output carries the structure it was run on, so the builder compares it against the displaced supercell of its position and names the directory and the displacement on a mismatch.

Zero-padded names (grid-001, disp-001) get both orderings right in the first place, whether the shell expands a glob lexicographically or a brace range in order. Scripts 1 and 2 write grid-{idx:03d} and disp-{k:03d}, and phonopy -d does the same. Padding is a convenience rather than a requirement: a wrong order is reported, not turned into force constants from mismatched forces.

The examples write the range out as grid-{001..025} rather than grid-*. The count is then visible in the command, and a missing grid point stops the builder on a path that does not exist instead of quietly passing a shorter list. grid-* also works, and is the shorter thing to type once the grid is known to be complete.

What the file holds#

aniso_qha_dataset.hdf5 is self-contained. Per grid point it stores the relaxed cell, the supercell and primitive matrices, the raw displacements and forces, the static internal energy $U$, and optionally the electronic states. The displacements and forces are kept in phonopy’s native displacement-force dataset form – type-1 (one displaced atom per supercell) or type-2 (dense/random) – and tagged with which one, so the reader picks the force-constant solver from the type instead of inferring it. Storing them raw rather than as force constants keeps the file independent of the force-constant method, and makes it an archive that outlives the calculator scratch.

The electronic states always carry the eigenvalues, the k-point weights and the electron count. They carry the k points and the mesh as well, but only when the static calculation sampled a Gamma-centred regular mesh. Those last two are what the linear tetrahedron method needs, so a file without them is integrated only by the k-point sum. The two readers then part company: the run_anisotropic_qha API falls back on the k-point sum, while the phonopy-anisotropic-qha command stops instead. Step 4 gives the reason.

The layout, as written by phonopy-anisotropic-qha-dataset and read back by phonopy.qha.anisotropic_dataset.read_aniso_qha_dataset:

/                    attrs: creator, phonopy_version, calculator, length_unit,
                            crystal_system, free_dof, tie_description,
                            n_grid_points, grid_shape
/grid/NNN            attrs: index, internal_energy, displacement_type
    lattice                          (3, 3)       the relaxed grid-point cell
    lattice_lengths                  (3,)         a, b, c
    scaled_positions                 (natom, 3)
    numbers                          (natom,)
    masses                           (natom,)
    magnetic_moments                 (natom,)     only when the cell has them
    supercell_matrix                 (3, 3)
    primitive_matrix                 (3, 3)
    displaced_atoms                  (ndisp,)     type-1 only
    displacements                    (ndisp, 3)   type-1
                                     (ndisp, natom_super, 3)   type-2
    forces                           (ndisp, natom_super, 3)
    electronic_states/eigenvalues    (nspin, nkpt, nband)   optional
    electronic_states/weights        (nkpt,)
    electronic_states/n_electrons    scalar
    electronic_states/fermi_energy   scalar
    electronic_states/spin_degeneracy  scalar     only when set
    electronic_states/kpoints        (nkpt, 3)    only with a regular mesh
    electronic_states/mesh           (3,)         only with a regular mesh

grid_shape is present only when the cells reached the builder as a tensor grid; --compare-eos needs it. NNN is the position in the --static list, and the analysis reads the lattice parameters from lattice rather than from that number.

4. Run the anisotropic QHA#

Run the analysis directly on the intermediate dataset:

% phonopy-anisotropic-qha aniso_qha_dataset.hdf5 --tmax 1000 --dt 10 \
    --contour-temp 0 500 1000 --compare-eos

--tmax and --dt are in K and set the temperature grid, here 0 to 1000 K in steps of 10 K. The thermal expansions are central differences on that grid, so the top temperature is consumed and the results stop one step below --tmax: 990 K for the command above.

--mesh sets the phonon sampling mesh and defaults to 200, denser than run_qha’s 100. The axial thermal expansions need the denser mesh, while the volumetric expansion is already converged at 100.

The command excludes the three acoustic modes at $\Gamma$ from the phonon thermal properties. --no-exclude-gamma-acoustic includes the ones whose frequencies are positive, as phonopy did before this option existed. See EXCLUDE_GAMMA_ACOUSTIC.

The command rebuilds one Phonopy per grid point, with force constants from the stored displacements and forces, runs run_anisotropic_qha, and writes lattice_parameters-temperature.dat, axial_thermal_expansion.dat, volume-temperature.dat and anisotropic_qha.png.

With exactly two free lattice DOF it also writes the F(a, c) contour maps. --decompose-contours adds the $U$ / $F_\mathrm{ph}$ / $F_\mathrm{el}$ / total panels, which say which term makes the valley and which one moves it: $U$ carries almost the whole curvature, while $F_\mathrm{ph}$ and $F_\mathrm{el}$ are nearly flat ramps whose tilt is what walks the minimum as the temperature rises.

The contour axes are the strain from the lattice that minimizes $U$ rather than the lattice parameters themselves, so that two calculations can be put side by side. That origin is a property of the electronic-structure calculation alone, while the lattice located at $T = 0$ already carries the zero-point term and moves when the vibrational model does. The window is half again as wide as the sampled cells, --margin, so that the located minimum stays visible when it sits near the edge. Its size is fixed by the grid and not by where the reference falls, so two calculations drawn with the same margin come out at the same scale; a reference far from the centre of the grid can push some sampled cells out of view. The levels are chosen geometrically from the range of the data and rounded, so they read as 0.2, 0.5, 1, 2, 5 rather than as whatever the maximum happens to be divided into forty. --plot-format pdf writes the figures as PDF instead of PNG.

Each contour figure is written with a .dat beside it holding the fit it draws: the polynomial coefficients, the centre and scale they are defined against, and the sampled cells with the value fitted at each. A contour map is a picture of ten numbers, and those ten numbers replot at any resolution and in any style, which a bitmap does not.

--compare-eos adds a volume-path cross-check along the main diagonal of the grid. The diagonal comes from the grid shape the builder recorded, so the cross-check needs a --grid run in step 1 whose cells were gathered in order. Without a recorded shape the command skips the cross-check and lists the cells by volume with their $c/a$, so that --eos-index can name a path by hand. Cells that share a $c/a$ differ in volume alone, and differing in volume alone is what a volume-path equation-of-state fit assumes. Any five or more such cells are printed as a ready-made --eos-index line.

The electronic free energy $F_\mathrm{el}$ is added whenever the dataset carries the electronic states, and --no-electronic leaves it out. The integration is the linear tetrahedron method. That method needs the k-point grid the states were computed on, and without that grid the command stops. It does not fall back on the k-point sum. The sum integrates a delta-function density of states, and so needs far more irreducible k points than a mesh chosen for the total energy has. The run_anisotropic_qha API is looser and does fall back on the sum, one grid point at a time, so a script calling the API has to check the mesh for itself. The run prints the energy window and the grid spacing it integrated over. The result is written to fel.hdf5, which --electronic-free-energies takes back on a later run, since the integration is the same every time.

Converge the mesh on the thermal expansion rather than on $F_\mathrm{el}$. Two meshes can agree closely on $F_\mathrm{el}$ and still differ severalfold in $\alpha_c$, because the expansion is a derivative of the free-energy surface and the error varies across the lattice grid.

The same analysis, driven from the API:

Script 3 – the analysis through run_anisotropic_qha#

import numpy as np
from phonopy import run_anisotropic_qha
from phonopy.qha import anisotropic_output, anisotropic_plot
from phonopy.qha.anisotropic_dataset import read_aniso_qha_dataset

dataset = read_aniso_qha_dataset("aniso_qha_dataset.hdf5")

phonopys = []
internal_energies = []
electronic_structures = []
for point in dataset.grid_points:
    # to_phonopy() rebuilds the Phonopy and force constants from the stored
    # dataset, picking the site-symmetry or symfc solver by dataset type.
    phonopys.append(point.to_phonopy())
    internal_energies.append(point.internal_energy)
    electronic_structures.append(point.electronic_states)

has_electronic = all(e is not None for e in electronic_structures)
temperatures = np.arange(0, 1001, 10.0)  # one extra point for finite diff
result = run_anisotropic_qha(
    phonopys,
    temperatures,
    internal_energies=internal_energies,
    electronic_structures=electronic_structures if has_electronic else None,
    # mesh defaults to 200 here, denser than run_qha's 100: the axial split
    # needs it, while the volumetric expansion is converged at 100.
)

anisotropic_output.write_lattice_parameters_temperature(result)
anisotropic_output.write_axial_thermal_expansion(result)
anisotropic_output.write_volume_temperature(result)

fig = anisotropic_plot.plot_anisotropic_qha(result)
fig.savefig("anisotropic_qha.png")

# The F(a, c) diagnostics, one file per temperature. plot_component_contours
# splits the surface into its U / F_ph / F_el / total parts, which is what
# shows where the valley comes from; both need exactly two free lattice DOF
# and return the names they wrote.
contour_temperatures = [0.0, 300.0, 600.0, 1000.0]
anisotropic_plot.plot_F_contours(result, contour_temperatures)
anisotropic_plot.plot_component_contours(
    result,
    internal_energies,
    electronic_structures if has_electronic else None,
    contour_temperatures,
)

Both forms run the same analysis. run_anisotropic_qha detects the free lattice DOF from the input cells, and then, at each temperature:

Add the terms at every sampled cell, $F_i(T) = U(a_i, c_i) + F_\mathrm{ph}(a_i, c_i; T)$, with $F_\mathrm{el}(a_i, c_i; T)$ when it is included and a $pV_i$ term when a pressure is given.
Fit a polynomial of total degree $n$ in the free lattice parameters to those $F_i(T)$, by least squares. --polynomial-degree sets $n$, 3 by default. For $d$ free DOF the polynomial has $\binom{n + d}{n}$ terms – 10 for two free DOF at degree 3 – and the grid needs at least that many cells, or the fit is rank deficient and says so.
Minimize that polynomial, from the centroid of the sampled cells and using its analytic gradient. The minimizer is $a(T)$, $c(T)$. A minimum outside the sampled range is reported as an extrapolation.

The axial thermal expansions then follow by central differences over the temperature grid,

\[\alpha_a = \frac{1}{a}\frac{da}{dT}, \qquad \alpha_c = \frac{1}{c}\frac{dc}{dT},\]

The fit is over the lattice only: each temperature is minimized on its own, and the raw $a(T)$, $c(T)$ are differenced. The calculator’s harmonic free energies carry no sampling scatter, so differencing them directly is the right thing to do here. A free energy that does carry scatter needs the lattice parameters smoothed along temperature first, which is --smooth-lattice; step 5 covers it, since that is where such free energies come from.

The energies are normalized per primitive cell, in eV: the internal energies, the phonon free energies and the electronic free energies. The volumes the analysis works in are that cell’s too, which is what a pressure run builds its $pV$ term on. Which cell it is comes from primitive_matrix in phonopy_disp.yaml, which phonopy-init sets to auto unless --pa says otherwise.

Nothing else on the page follows the primitive cell. The grid, the lattice parameters and the axial thermal expansions are all the conventional cell’s, as step 0 says.

The static single point need not be run on that cell. Either the conventional unit cell or the primitive cell will do, whichever suits the calculation – a long rhombohedral cell is easier to handle in its hexagonal setting, while a large conventional cell may be cheaper to run as its primitive cell. The builder settles the normalization: it reads the number of atoms in the cell $U$ was computed on, and multiplies $U$ by the primitive cell’s share of it. A static point already on the primitive cell is left alone. One run on a cell holding four of them is divided by four, and the builder says so:

The static single point was run on a cell holding 4 primitive cells, so U is
stored as 1/4 of the calculator energy, matching the phonon free energy.

Which cell it was is read off the atom counts rather than assumed, because strain changes the volumes but not them.

Two restrictions come with this freedom. The static cell has to be the conventional unit cell when --phonon is omitted, since that cell then becomes the grid point and the free lattice DOF are read from its rows; the builder stops and says so. And with --phonon given, the static cell must be a cell of the grid point it is paired with, either of the two. A cell from another grid point matches neither, so a mis-pairing is still caught.

$F_\mathrm{el}$ is normalized the same way, where the analysis integrates it from the stored electronic states. The states record the cell they were computed on, so states already on the primitive cell are not scaled a second time. Arrays handed to run_anisotropic_qha directly are the exception: internal_energies, electronic_free_energies and phonon_free_energies are taken as they are, and have to be per primitive cell already.

Supplying $F_\mathrm{el}$ ready-made#

The tetrahedron integration costs of order a minute per grid point on a dense mesh, and the analysis works through the grid points one after another. Computing $F_\mathrm{el}$ outside lets one result be reused across runs, and lets the grid points be spread over processes or jobs. Pass it as electronic_free_energies, the counterpart of phonon_free_energies for the electronic term.

Script 4 – $F_\mathrm{el}$ computed apart from the analysis#

import numpy as np

from phonopy import run_anisotropic_qha
from phonopy.qha.anisotropic_dataset import read_aniso_qha_dataset
from phonopy.qha.electron import compute_free_energy_by_tetrahedron

dataset = read_aniso_qha_dataset("aniso_qha_dataset.hdf5")
temperatures = np.arange(0, 1001, 10.0)  # one extra point for finite diff

fe_el = np.column_stack(
    [
        # F_el(T) - F_el(0) in eV per primitive cell, one grid point at a time.
        compute_free_energy_by_tetrahedron(point.electronic_states, temperatures)[0]
        for point in dataset.grid_points
    ]
)

result = run_anisotropic_qha(
    [point.to_phonopy() for point in dataset.grid_points],
    temperatures,
    internal_energies=[point.internal_energy for point in dataset.grid_points],
    electronic_free_energies=fe_el,
)

The grid points are independent, so this loop splits across processes or jobs unchanged; each one writes its own column of fe_el. The integration itself is already threaded through BLAS, so cap the threads per process when several run on one node.

The values are anchored at T = 0 and normalized per primitive cell, consistently with internal_energies. electronic_structures and electronic_free_energies are two ways of giving the same term, so pass one or the other.

phonopy-anisotropic-qha takes the same thing from a file. Write it as an ElectronicFreeEnergies, which carries the temperatures it was computed on:

from phonopy.qha.free_energy_io import (
    ElectronicFreeEnergies,
    write_free_energies_hdf5,
)

write_free_energies_hdf5(
    ElectronicFreeEnergies(
        temperatures=temperatures,
        free_energies=fe_el,
        # Optional, and what lets the command check the file against the grid
        # it is used with.
        lattice_lengths=np.array(
            [np.linalg.norm(point.cell.cell, axis=1) for point in dataset.grid_points]
        ),
    ),
    "fel.hdf5",
)

write_free_energies_hdf5 writes the phonon terms as well, and the file records which term it holds, so reading it back as another one is refused rather than silent.

% phonopy-anisotropic-qha aniso_qha_dataset.hdf5 --tmax 1000 --dt 10 \
    --electronic-free-energies fel.hdf5

The file carries the temperatures it was computed on. The command compares them with the grid --tmax and --dt ask for and stops if they differ, instead of pairing the rows as they come. --electronic-free-energies replaces --electronic; passing both stops the command. It also skips --compare-eos, since the volume-path driver takes the electronic term as states and the two paths would then carry different physics.

5. Variant: temperature-dependent force constants from per-grid-point MLPs#

The displacements of step 2 are small, and the forces they give are almost harmonic. Anharmonic effects appear only in supercells displaced further than that. Here such supercells are made by moving every atom at once, by amounts drawn at random from the thermal distribution of the harmonic crystal at a chosen temperature.

The force constants of step 2 define that harmonic crystal. It separates into independent normal modes. Each mode is a harmonic oscillator, so its own normal coordinate is Gaussian about zero, and the width of that Gaussian is set by the mode’s frequency and by the temperature. One snapshot draws every mode from its own Gaussian, and the displacement of each atom follows from the sum over the modes. “The thermal distribution” below writes the width down.

The calculator then computes forces on supercells displaced that way at that temperature. These displacements are usually not small, so the forces carry more of the anharmonic contribution than those of step 2 do, and are enough to train an MLP on.

The displacements, the forces and the supercell energies together are used to build the MLP at every grid point. Each MLP is fitted to that grid point’s own training set.

The MLP is then used, within the temperature range its training set covers. An MLP is an intermediate representation, convenient because it interpolates between the temperatures it was trained at. The self-consistent harmonic approximation (SSCHA) runs with it at every grid point and every temperature, and returns force constants that change with temperature. The anharmonic free energies follow from those. The analysis takes them and minimizes $F(a, c; T)$ over the lattice at each temperature. What comes out is $a(T)$, $c(T)$ and the axial thermal expansions.

Where this step fits#

Steps 0 to 4 displace one atom at a time and need a few supercells per grid point. They give the quasi-harmonic answer.

Step 5 needs a training set at every grid point instead. In the script below that is 50 structures at each of four temperatures, so 200 calculator runs per grid point. What it gives back is force constants that change with temperature. Use it when anharmonicity is large enough to matter for the property being computed.

Steps 0 to 3 are as in the Overview, up to and including the dataset. What differs comes after it: the MLPs and the free energies,

        flowchart TD
    DS(["aniso_qha_dataset.hdf5<br/>cells, U, F_el,<br/>harmonic force constants"])
    DS --> DISP["thermal displacements<br/>one shared draw of<br/>standard normals"]
    DISP --> CALCT{{"calculator forces"}}
    CALCT --> TR(["merged.yaml<br/>per grid point"])
    TR --> DEV["train one MLP<br/>per grid point"]
    DEV --> MLP(["polymlp.yaml<br/>per grid point"])
    MLP --> SSCHA["SSCHA<br/>per grid point and temperature"]
    TR --> SSCHA
    SSCHA --> FE(["F_ph(T) per<br/>grid point"])

and the analysis they meet in:

        flowchart TD
    DS2(["aniso_qha_dataset.hdf5"])
    FE2(["F_ph(T) per<br/>grid point"])
    DS2 --> AQ["run_anisotropic_qha"]
    FE2 -->|"phonon_free_energies"| AQ
    AQ --> RES(["a(T), c(T),<br/>alpha_a, alpha_c"])

One MLP per grid point#

The MLP gives the temperature-dependent force constants. The SSCHA step below computes the anharmonic phonon free energy $F_\mathrm{ph}(a, c; T)$ from them.

$U(a, c)$ comes from the calculator on the static grid of step 1, not from the MLP. The analysis minimizes $F$ over the lattice at each temperature. The shape of $U$ near the minimum therefore decides where the minimum lies. The axial expansions measure how that minimum moves as the temperature rises. A small error in the shape of $U$ moves the minimum, and the axial expansions follow that error. That is why $U$ is taken from the calculator rather than from the MLP. The same sensitivity is behind the choice of the linear tetrahedron method for $F_\mathrm{el}$ in step 4: both terms shape the surface that the analysis differentiates.

Each MLP is evaluated only at its own grid point. It has to reproduce the forces and energies of displaced supercells there and nowhere else. The forces give the temperature-dependent force constants. The energies enter the SSCHA free energy as differences from the undisplaced supercell of that grid point. The grid as a whole carries the lattice dependence of $F_\mathrm{ph}$, so in this procedure the MLPs are never asked to interpolate in the lattice parameters.

The MLPs are fitted independently, so their errors can jump from one grid point to the next. Drawing the displacements from one set of standard normals is expected to smooth those errors out. See Sharing one draw of standard normals.

The thermal distribution#

The training structures are supercells displaced along the thermal distribution of the harmonic crystal. Drawing them costs nothing beyond the force constants of step 2, and the amplitudes are the ones the crystal actually visits at the temperature asked for.

The distribution is drawn once, from the harmonic force constants. SSCHA iterates its own distribution to self-consistency with the anharmonic force constants; this draw does not. In practice that has been enough for a training set.

The thermal distribution writes down the variance $\sigma_{\mathbf{q}\nu}^2$ of each mode and the displacements that follow from one draw of standard normals $\xi_{\mathbf{q}\nu}$. A grid point’s own frequencies and eigenvectors enter there, which is what Sharing one draw of standard normals below rests on.

The training displacements#

Four steps, once the dataset of step 3 exists. The harmonic force constants it carries are what set the widths of the draw.

Pick temperatures covering the temperatures script 7 will run. An MLP tends to be poor outside the range it was trained on. Script 5 below trains at 0, 100, 250 and 400 K, and script 7 runs up to 400 K, so every SSCHA run stays inside the trained range.
Generate the displacements with the script below. It reads the harmonic force constants from the dataset of step 3, and uses one draw of standard normals at every grid point; see Sharing one draw of standard normals.
Run the calculator on every supercell, and collect the forces of each set with phonopy-init -f. This writes one phonopy_params.yaml per set, holding its displacements, forces and supercell energies.
Merge the phonopy_params.yaml files of each grid point, one per temperature, into a single training set, and train that grid point’s MLP on it.

Script 5 – the thermal training displacements#

from pathlib import Path

import numpy as np

from phonopy.interface.vasp import write_vasp
from phonopy.qha.anisotropic_dataset import read_aniso_qha_dataset

TRAIN = Path("train")
TEMPERATURES = (0.0, 100.0, 250.0, 400.0)  # K
SNAPSHOTS = 50  # structures per grid point and temperature
SEED = 20260815

dataset = read_aniso_qha_dataset("aniso_qha_dataset.hdf5")
TRAIN.mkdir(parents=True, exist_ok=True)

for temperature in TEMPERATURES:
    normals = None
    for point in dataset.grid_points:
        phonon = point.to_phonopy()
        phonon.init_random_displacements()
        rd = phonon.random_displacements
        if normals is None:
            # Drawn once and reused at every grid point; the shapes are set by
            # the supercell matrix and the primitive cell, which they share.
            normals = rd.draw_standard_normals(
                SNAPSHOTS, random_seed=SEED + int(temperature)
            )
            # Keep the draw itself. The displacements below record the
            # ensemble, but only these reproduce it, and extending it with
            # first_snapshot needs them.
            np.savez_compressed(
                TRAIN / f"normals-{int(temperature)}K.npz",
                ii=normals[0],
                ij=normals[1],
                seed=SEED + int(temperature),
                temperature=temperature,
            )
        rd.run(temperature, standard_normals=normals)
        phonon.dataset = {"displacements": rd.u.copy()}

        # The stored index is 0-origin; the directories are numbered from 001.
        set_dir = TRAIN / f"grid-{point.index + 1:03d}-{int(temperature)}K"
        set_dir.mkdir(parents=True, exist_ok=True)
        # The displacements have to be saved beside the supercells: it is what
        # phonopy-init -f attaches the forces to below.
        phonon.save(
            set_dir / "phonopy_disp.yaml",
            settings={"force_constants": False, "displacements": True},
        )
        for i, cell in enumerate(phonon.supercells_with_displacements, 1):
            disp_dir = set_dir / f"disp-{i:03d}"
            disp_dir.mkdir(exist_ok=True)
            write_vasp(disp_dir / "POSCAR", cell)

Each set directory then holds phonopy_disp.yaml and disp-001/POSCAR .. disp-050/POSCAR, the same layout as the phonon grid of step 2, and train/ holds one normals-*.npz per temperature beside them. Run the calculator in every disp-*, then collect the forces of each set:

% phonopy-init -f disp-*/vaspout.h5 --save-params
# -> phonopy_params.yaml, with the displacements, forces and supercell energies

A grid point has one such set per temperature, and its MLP is trained on all of them at once. The sets are merged by interleaving, so that the temperatures alternate through the merged list instead of following one another in blocks.

Merging the temperatures gives the reason, which is how ntrain and ntest cut the merged list:

Script 6 – one training set per grid point, from its temperatures#

from pathlib import Path

import numpy as np

import phonopy

TRAIN = Path("train")
TEMPERATURES = (0, 100, 250, 400)
N_GRID = 25

for index in range(1, N_GRID + 1):
    sets = [
        phonopy.load(
            TRAIN / f"grid-{index:03d}-{t}K" / "phonopy_params.yaml",
            produce_fc=False,
            log_level=0,
        )
        for t in TEMPERATURES
    ]

    # Interleave the temperatures, so that any prefix and any suffix of the
    # merged set holds them in equal parts.
    merged = {}
    for key in ("displacements", "forces", "supercell_energies"):
        stacked = np.array([s.dataset[key] for s in sets])
        merged[key] = stacked.swapaxes(0, 1).reshape(-1, *stacked.shape[2:])

    mlp_dir = TRAIN / f"grid-{index:03d}"
    mlp_dir.mkdir(parents=True, exist_ok=True)
    phonon = sets[0]
    phonon.dataset = merged
    phonon.save(
        mlp_dir / "merged.yaml",
        settings={"force_sets": True, "displacements": True},
    )

To check the sets before starting the calculator, see checking the training displacements.

Then train one MLP per grid point on its merged set. phonopy --pypolymlp always writes the MLP as polymlp.yaml in the current directory, and the name cannot be changed from the command line. So run the training inside that grid point’s own directory. Run it in train/ instead and the 25 grid points write over one file, leaving only the last one:

% cd train/grid-013
% phonopy merged.yaml --pypolymlp --mlp-params="ntrain=..., ntest=..." -v
# -> train/grid-013/polymlp.yaml

Script 7 reads the MLPs back from train/grid-NNN/polymlp.yaml, which is where this puts them, and the cell and the starting force constants from the merged.yaml beside them.

The displacements are drawn at random, so a training set of a given size is one draw among many. A second draw of the same size would give a different training set, and with it a different MLP and different results downstream.

That difference can be measured. Pick one grid point. Draw a second set of the same size there with a different seed, and train a second MLP on it. Comparing the phonon frequencies the two MLPs give is the cheap check, but it may not be enough. The frequencies can agree closely while the axial expansions still differ, so compare the quantity you intend to report.

What the draw leaves at its defaults#

Script 5 leaves the other parameters of init_random_displacements at their defaults, and what the draw leaves at its defaults says which of them change the displacements.

A quasi-harmonic grid can reach grid points that are dynamically unstable, and the draw takes $|\omega|$, so an imaginary mode there is drawn as a real mode of the same magnitude. Look at the frequencies of step 2 before training on such a grid point.

The descriptor and the amount of training data#

The descriptor and the amount of training data covers what --mlp-params sets, what a descriptor costs to evaluate, and how the ridge penalty pypolymlp selects says whether the training set is large enough for the descriptor.

One ladder of descriptors lists the feature counts and relative evaluation times of nine descriptors for a one-element system.

The SSCHA evaluates the descriptor once per snapshot per iteration, and here that cost is paid at every grid point and every temperature. Choose the descriptor and the training-set size from fits made at one grid point, before the sweep is started.

A lower force RMSE does not have to reach the thermal expansion, and the thermal expansion is what this step is for. “Validate the MLP” below says what to compare instead.

Pinning the ridge penalty across the grid#

pypolymlp fits every penalty in reg_alpha_params and keeps the one with the smallest test RMSE, chosen separately at every grid point. The analysis differentiates across the grid, so a penalty that changes from one grid point to the next puts a step into the quantity being differentiated. Collapse the range to one point to pin it:

% phonopy grid-NNN-merged.yaml --pypolymlp \
    --mlp-params="ntrain=..., ntest=..., reg_alpha_params = -3.0 -3.0 1" -v

The three numbers are linspace(p0, p1, p2) of the base-10 logarithm, so -3.0 -3.0 1 is alpha = 1e-3 alone, against the default -3.0 1.0 5 of 1e-3 to 1e1 in five steps.

Validate the MLP#

Validating the MLP compares the MLP forces against the calculator forces on the structures held out as the test set, one temperature at a time. The thermal supercells of this step already carry calculator forces, so nothing new has to be run with the calculator.

Make that comparison at a few grid points. The MLPs are fitted independently, so an MLP that is weak at one grid point says nothing about the others.

Computing the free energies with SSCHA#

$F_\mathrm{ph}$ is the SSCHA free energy defined in SSCHA, computed from force constants that change with temperature. It is no longer the harmonic expression of “The free energy” above.

phonopy-mlpsscha computes the free energy at one (temperature, grid point) pair, from that grid point’s own training set and MLP. Script 7 below runs the temperatures of one grid point, and script 8, listed in the appendix at the end of this page, gathers what the runs wrote. Save them as script7.sh and script8.py. The axial thermal expansions are computed from the free energies of all the pairs.

The temperatures must be decided before the first SSCHA run, since the gather places every run on one temperature grid. The two scripts carry 0 to 400 K in 10 K steps. These 41 temperatures over the 5 x 5 lattice grid of step 1 are 1,025 runs.

$a(T)$ and $c(T)$ are obtained by minimizing the free-energy surface at each temperature. The axial expansions then come from an Einstein fit of $a(T)$ and $c(T)$ over the whole temperature range, described in running the analysis. A temperature at the end of that range is not interpolated, so it is recommended to run above the highest temperature to report.

One SSCHA run reads train/grid-NNN/merged.yaml for the cell and train/grid-NNN/polymlp.yaml for the forces, and writes what every iteration sampled to its own file. Script 7 is one grid point at every temperature, which is the unit a job usually gets:

Script 7 – the SSCHA runs of one grid point#

#!/bin/bash
# One grid point at every temperature:  ./script7.sh 13
set -e
g=$(printf '%03d' "$1")  # 13 -> 013

for t in {0..400..10}; do
    phonopy-mlpsscha "train/grid-$g/merged.yaml" \
        --mlp "train/grid-$g/polymlp.yaml" \
        -t "$t" \
        --snapshots 2000 \
        --iterations 16 \
        --mesh 200 \
        --random-seed 1000 \
        -o "sscha-g$g-t${t}K.hdf5"
done

The argument is the grid point, numbered from 1 as the directories are. A run reads those two files and nothing else: aniso_qha_dataset.hdf5 is read only by the gather, for the lattice lengths it places the runs by.

The force constants a run starts from are those of merged.yaml, which phonopy fits with symfc to the training displacements of all four temperatures. They are not the harmonic force constants of step 2: the training set was drawn at 0, 100, 250 and 400 K, so the fit already carries some of the thermal displacement. How far they sit from the self-consistent force constants of a given temperature is what choosing the transient reads off the listing.

Gathering the runs below covers how the calls are spread over jobs and how script 8 puts them back together. It is worth making one of them by hand first, with -v added, which lists its iterations. The transient is chosen from that listing, and -vv adds the force-constant fit.

Sampling and averaging are separate steps. The averaging is where the transient is chosen, and doing it apart means choosing again costs a second of arithmetic rather than the whole sweep.

The options of an SSCHA run and its error#

The options of a run covers what each option sets, and the error of the average covers how the error of the mean follows from --snapshots and the number of iterations kept. Script 7 draws 2000 supercells against phonopy’s own 1000, and makes 16 iterations against its 10.

The --mesh of script 7 matches the --mesh of the analysis, which keeps one sampling through the calculation.

A run stores every iteration and averages none of them. Script 8 takes the mean over the iterations after the transient, one of them by default, and --transient sets how many. Choosing another is a second gather and costs no sampling, and --transient on a run marks its listing alone.

Gathering the runs into `fph.hdf5`#

One SSCHA run is minutes with the lightest descriptor and longer with a heavy one, and the 1,025 of them are independent of each other, so they can be split over processes, nodes or jobs however is convenient. Spreading the sweep over a cluster is one way. The same --random-seed draws the same supercells, so a run does the same sampling wherever it is made and however often it is repeated; what a rerun reproduces says where that stops.

% ./script7.sh 13
Wrote sscha-g013-t0K.hdf5
Wrote sscha-g013-t10K.hdf5
...

The gather places each run at the one of its TEMPERATURES the run’s own temperature matches, to within 1e-3 K. A submitting script can therefore carry its own temperatures, since one that is off by rounding still lands on the temperature the gather expects. A run at a temperature that is on no such grid point is left out, so a sweep computed over a wider or a finer grid than TEMPERATURES gathers to TEMPERATURES.

Sampling writes one sscha-g*K.hdf5 per run. Averaging writes fph.hdf5, which is the file the analysis reads.

The two are different types, not two spellings of one. A sscha-g*K.hdf5 file holds an SSCHATrace, every iteration and no average, and is what another transient is taken from. fph.hdf5 holds SSCHAFreeEnergies, the averages and the transient they were taken with. Handing the analysis a run stops with a message about its type, rather than averaging it over iterations nobody chose.

Script 8 gathers the sscha-g*K.hdf5 files into fph.hdf5:

% python script8.py
Wrote fph.hdf5, 41 temperature(s) x 25 grid point(s), from 1025 file(s)

Each file is placed by the lattice lengths and the temperature it carries rather than by its name, so the order they are gathered in does not matter.

Choosing the transient#

--transient sets how many iterations are left out of the averages. The default is DEFAULT_TRANSIENT, 1. Iteration 1 uses the force constants fitted to merged.yaml, so its free energy is that of those force constants and not of self-consistent ones.

Reading the run says how the listing -v prints is read and when to raise --transient.

Check the lowest temperature as well as the highest. The starting force constants are one fit over the whole training range, so they sit furthest from the self-consistent ones at its two ends.

The run files hold every iteration and no average. Script 8 computes the averages, leaving out the first --transient iterations of each run:

% python script8.py --transient 2

Trying --transient 3 will simply overwrite the existing fph.hdf5.

The runs gathered into one file must all have made the same number of iterations. Averaging over different numbers would weight the grid unevenly, so the gather stops and names the numbers it found.

Supplying the free energies to the analysis#

In steps 0 to 4 each Phonopy instance of phonopys carries one force-constant array, and it serves every temperature. The analysis computes $F_\mathrm{ph}$ from them:

run_anisotropic_qha(phonopys, temperatures, internal_energies=internal_energies)

Here they are computed elsewhere and passed to run_anisotropic_qha instead:

run_anisotropic_qha(
    phonopys,
    temperatures,
    internal_energies=internal_energies,
    phonon_free_energies=free_energies,  # shape (temperatures, grid points)
)

The Phonopy instances then supply only the cells and volumes. They can be built without force constants, which are never read. Without phonon_free_energies, run_anisotropic_qha computes the phonon free energy itself, and the force constants have to be there.

The values are per primitive cell. They are normalized the same way as internal_energies, and they do not include the static energy, which internal_energies already carries.

An MLP evaluates the undisplaced supercell, and its free energies are measured from that energy. SSCHAFreeEnergies records it as reference_energies, and phonopy-anisotropic-qha --use-mlp-internal-energies adds it back and takes U = 0, which puts the whole surface on the potential’s own energy scale rather than the calculator’s.

Since no force constants are read, the analysis also runs on a dataset built from the static grid alone. That is all a method has to work with when it never computed calculator phonons:

% phonopy-anisotropic-qha-dataset phonopy_disp.yaml \
    --static static-grid/grid-{001..025}/ -o aniso_qha_dataset.hdf5
No --phonon given: building from the static grid alone, 25 grid point(s), with
no displacements or forces. Such a dataset is for use with the
phonon_free_energies argument of run_anisotropic_qha.
  grid 1 U=... eV n_disp=0
  ...
  grid 25 U=... eV n_disp=0

Its grid points carry the cells, $U$ and the electronic states, but no harmonic force constants. The SSCHA runs of step 5 read merged.yaml and the gather reads only the cells, so both are content with such a dataset. What it cannot do is the draw of step 5: script 5 takes the widths of the thermal distribution from those harmonic force constants, and there are none here.

Running the analysis#

The file is given to the analysis with --phonon-free-energies:

% phonopy-anisotropic-qha aniso_qha_dataset.hdf5 --phonon-free-energies fph.hdf5

The file carries the temperatures it was computed on, and the command takes its temperature grid from them unless --tmax or --dt says otherwise. It also carries the lattice lengths of the grid points, when they were written, and those are checked against the dataset. A file computed on another machine cannot then be paired with the wrong grid.

A free energy from a sampled method carries the scatter of its sampling. The lattice parameters that minimize the free-energy surface at one temperature move with that scatter, and they move independently of those at the next temperature. Without smoothing the analysis takes each of $da/dT$, $db/dT$ and $dc/dT$ as a central difference between neighboring temperatures, which amplifies that scatter. The command above therefore smooths the lattice parameters before differentiating them, and says so as it runs.

The form fitted to each lattice parameter is

\[a(T) = a_0 + \sum_i A_i \frac{\theta_i}{e^{\theta_i / T} - 1},\]

with amplitudes $A_i$ of opposite sign and Einstein temperatures $\theta_i$, and the same for $b(T)$ and $c(T)$. Each term is zero at $T = 0$ with zero slope, so the fitted expansion vanishes at 0 K, as the third law requires. A general-purpose smoother assumes only smoothness, and a spline fitted to the same data returns a finite expansion at 0 K instead. Other forms share the property; this one is also the usual model for a lattice parameter that contracts at low temperature and expands at high temperature.

Differentiating it term by term gives

\[\frac{da}{dT} = \sum_i A_i \frac{(\theta_i / 2T)^2}{\sinh^2(\theta_i / 2T)},\]

which is the derivative the axial expansions are built from once the fit is made.

The smoothing is --smooth-lattice einstein, the default whenever --phonon-free-energies is given. --smooth-terms sets how many terms, 2 by default. More terms follow a curve more closely, and follow its scatter more closely too. --smooth-lattice none turns the smoothing off.

The free energies of steps 0 to 4 carry no sampling scatter, so the default there is none and the central differences are used. --smooth-lattice einstein works there as well, for the analytic derivative in place of them.

The same result comes from the API, with phonon_free_energies and lattice_smoothing:

result = run_anisotropic_qha(
    phonopys,
    temperatures,
    internal_energies=internal_energies,
    phonon_free_energies=free_energies,
    lattice_smoothing="einstein",
)

With phonon_free_energies given, run_anisotropic_qha skips the mesh sampling and ignores mesh.

Fitting a sum of Einstein terms is not a linear least squares, and different starting values converge to different curves. A converged curve can be wrong in shape: monotone where the data contracts, or dipping several times deeper than the data does. The fit therefore starts from a set of starting values, drops the curves whose shape disagrees with the data in those ways, and keeps the closest of what is left. If nothing is left, the command stops rather than returning a curve of the wrong shape.

A fit is worth seeing against what it was fitted to. The result therefore keeps the surface minima in unsmoothed_lattice_parameters, whether or not they were smoothed, and a smoothed run writes them as three columns more in lattice_parameters-temperature.dat: temperature, the smoothed $a$, $b$, $c$, then the same three before the smoothing. With --smooth-lattice none the two are the same numbers.

It also writes lattice_smoothing.png, one column per free lattice DOF. The upper row is the fit as a line over the minima as dots, and the lower row is what is left over, fit minus minimum, in $10^{-4}$ angstrom. The residuals are what the fit is judged on. Sampling scatter shows as a band around zero as wide as the scatter is, while a fit of the wrong shape shows as an excursion over a range of temperature, and that is the case to raise --smooth-terms for.

An unsmoothed run has nothing to compare against, so --smooth-lattice none writes the four columns alone and no lattice_smoothing.png.

A smoothed result carries the fitted model itself, so it answers at temperatures the run did not visit. lattice_parameters_at, axial_thermal_expansions_at, thermal_expansion_at and equilibrium_volumes_at take one temperature or an array of them:

Evaluating the smoothed lattice between the temperatures of the run#

a, b, c = result.lattice_parameters_at([293.15])[0]
alpha = result.axial_thermal_expansions_at([100.0, 200.0, 300.0])

The temperatures must lie within the range the run covered. A sum of Einstein terms models the data it was fitted to and nothing beyond it: a term of very high Einstein temperature is flat over the fitted range and turns on above it, so the model interpolates and refuses to extrapolate. A temperature outside the range raises ValueError, and a result from --smooth-lattice none has no model to evaluate and raises as well.

Appendix: one ladder of descriptors#

The descriptor and the amount of training data gives the method and a table for a two-element system. The table here is one example for a one-element system, measured with pypolymlp 0.20.5. The last column is the time to evaluate the descriptor once, relative to the first row, measured in one execution.

features	model parameters added to `--mlp-params`	relative time
781	nothing; phonopy’s defaults	1.0
1,176	`gaussian_params2 = 0 7 15`	1.2
2,600	`gaussian_params2 = 0 7 15, gtinv_maxl = 12 12`	5.6
3,848	`gaussian_params2 = 0 7 15, gtinv_order = 4, gtinv_maxl = 16 12 4`	7.1
6,820	`model_type = 4`	1.4
13,920	`model_type = 4, gaussian_params2 = 0 7 15`	1.4
22,495	`model_type = 4, gtinv_order = 6, gtinv_maxl = 16 12 4 1 1`	6.4
27,664	`model_type = 4, gaussian_params2 = 0 7 15, gtinv_maxl = 12 12`	5.9
45,680	`model_type = 4, gaussian_params2 = 0 7 15, gtinv_order = 6, gtinv_maxl = 16 12 4 1 1`	7.7

Phonopy’s defaults are model_type = 3, max_p = 2, gtinv_order = 3, gtinv_maxl = 8 8, gaussian_params2 = 0 7 10 and cutoff = 8.0. Phonopy passes them to pypolymlp itself. The other rows change one or two of them.

The evaluation time follows gtinv_maxl rather than the feature count. Two pairs in the table differ in gtinv_maxl alone, 1,176 against 2,600 and 13,920 against 27,664, and both cost about four times more at 12 12 than at 8 8. Raising model_type or the number of gaussians multiplies the feature count instead, and adds a few tens of per cent to the time: 781 to 13,920 is eighteen times the features for 1.4 times the time.

The SSCHA of step 5 evaluates the descriptor once per snapshot per iteration, at every grid point and every temperature. Read the last column as what a descriptor costs there, and pick the descriptor from fits made at one grid point before the sweep is started.

Appendix: the gather script#

What gathering the runs calls. It reads the dataset of step 3 for the lattice lengths it places the runs by, and nothing else of it. --transient says how many iterations at the start of each run to leave out of the averages.

Script 8 – the runs gathered into the free energies of the analysis#

import argparse
import glob

import numpy as np
from phonopy.qha.anisotropic_dataset import read_aniso_qha_dataset
from phonopy.qha.free_energy_io import (
    assemble_sscha_free_energies,
    write_free_energies_hdf5,
)
from phonopy.sscha.trace import read_sscha_trace_hdf5

DATASET = "aniso_qha_dataset.hdf5"
TEMPERATURES = np.arange(0, 410, 10.0)  # 0 to 400 K in 10 K steps
DEFAULT_TRANSIENT = 1


def assemble(
    dataset,
    transient=DEFAULT_TRANSIENT,
    pattern="sscha-g*K.hdf5",
    filename="fph.hdf5",
):
    """Gather what the runs wrote into the free energies the analysis reads."""
    paths = sorted(glob.glob(pattern))
    if not paths:
        raise SystemExit(f"No file matches {pattern}.")
    points = dataset.grid_points
    try:
        free_energies = assemble_sscha_free_energies(
            [read_sscha_trace_hdf5(path) for path in paths],
            TEMPERATURES,
            np.array([np.linalg.norm(p.cell.cell, axis=1) for p in points]),
            transient,
        )
    except ValueError as error:
        raise SystemExit(str(error)) from error
    write_free_energies_hdf5(free_energies, filename)
    n_temperatures, n_points = free_energies.free_energies.shape
    print(
        f"Wrote {filename}, {n_temperatures} temperature(s) x "
        f"{n_points} grid point(s), from {len(paths)} file(s)",
        flush=True,
    )


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "--transient",
        type=int,
        default=DEFAULT_TRANSIENT,
        help="how many iterations at the start of a run are its transient "
        "and are left out of the averages (default: %(default)s)",
    )
    args = parser.parse_args()

    assemble(read_aniso_qha_dataset(DATASET), args.transient)


if __name__ == "__main__":
    main()

Appendix: checking the training displacements#

Script 5 writes one set of displaced supercells per grid point and temperature, into train/grid-NNN-TK/. Script 6 merges the sets of one grid point into train/grid-NNN/merged.yaml. The calculator is run on the POSCAR of each disp-*, and phonopy --pypolymlp trains on merged.yaml.

The checks in this appendix compare what script 5 and script 6 wrote against the distribution the supercells were drawn from. They take seconds. The checks on the displacements read phonopy_disp.yaml and the dataset of step 3, so run them before the calculator. Script 6 merges the forces with the displacements, so the check on the merged sets has something to read only after the calculator, and before then it reports every grid point as having no merged.yaml.

What each check compares#

One check compares the amplitude of a set against the temperature its directory is named after. Checking the training displacements writes down the ratio it prints, which reference_u2 and sample_u2 compute here from the force constants of the grid point the set belongs to. check_amplitudes prints the ratio for every set, and a set with the wrong label gives a ratio far from 1.

Wrong force constants are not what this catches. The draw and the reference use the same $\Phi$, so force constants that are wrong at a grid point change both by the same amount and the ratio still comes out 1. Catch those at step 2, from the frequencies.

The next check compares neighbouring grid points. Grid points $g$ and $g'$ draw from $\tilde{\rho}_{\Phi^{(g)}}(T)$ and $\tilde{\rho}_{\Phi^{(g')}}(T)$, whose force constants differ only by the strain between the two grid points. The grid points also share one draw of standard normals (Sharing one draw of standard normals), so at a fixed temperature and snapshot their displacements differ only through $\Phi^{(g)}$ and $\Phi^{(g')}$. correlation measures how close the two sets are as

\[\rho(g, g') = \frac{\sum_{n} \sum_{l\kappa j} u_{l\kappa j}^{(g,n)} u_{l\kappa j}^{(g',n)}} {\Bigl[ \sum_{n} \sum_{l\kappa j} \bigl( u_{l\kappa j}^{(g,n)} \bigr)^2 \sum_{n} \sum_{l\kappa j} \bigl( u_{l\kappa j}^{(g',n)} \bigr)^2 \Bigr]^{1/2}},\]

where $u_{l\kappa j}^{(g,n)}$ is snapshot $n$ of the set at grid point $g$, and every sum runs over the $N$ snapshots and the supercell. Shared normals put $\rho$ near 1, and the spacing of the grid sets how far below 1 it falls. check_neighbours prints $1 - \rho$.

Two draws that do not share their normals are uncorrelated, so $\rho$ is 0. A value of $1 - \rho$ near 1 means the normals were not shared.

read_set compares the basis vectors of each set against those of the grid point its directory is named after. The set carries its own unit cell in phonopy_disp.yaml, the cell script 5 built the supercells from, and the dataset of step 3 carries the relaxed cell of each grid point as point.cell. Both comparands are cell, the 3x3 matrix whose rows are the basis vectors, and they agree element by element to 1e-8 Angstrom when the directory holds the grid point it names. read_set also counts the snapshots of the set.

check_merged reads the merged set of script 6. It counts the structures and checks that the structures of each temperature are at the offsets k::len(TEMPERATURES) the merge writes them to. It closes with the number of merged sets it read, and names the grid points that have no merged.yaml.

The script#

Script 9 – checking the training displacements#

from pathlib import Path

import numpy as np

import phonopy
from phonopy.qha.anisotropic_dataset import read_aniso_qha_dataset

DATASET = "aniso_qha_dataset.hdf5"
TRAIN = Path("train")
TEMPERATURES = (0, 100, 250, 400)
SNAPSHOTS = 50


def reference_u2(point, temperatures):
    """Return <u^2> per temperature, from the force constants of one point.

    uu[i, j] is the 3x3 block of atoms i and j, in Angstrom^2.

    """
    phonon = point.to_phonopy()
    phonon.init_random_displacements()
    rd = phonon.random_displacements
    reference = {}
    for temperature in temperatures:
        rd.run_correlation_matrix(float(temperature))
        reference[temperature] = np.einsum("iiaa->", rd.uu)
    return reference


def sample_u2(displacements):
    """Return <u^2> over the snapshots of one set."""
    return (displacements**2).sum() / len(displacements)


def correlation(u_a, u_b):
    """Return the correlation of two sets of displacements."""
    a, b = u_a.ravel(), u_b.ravel()
    return float(a @ b / np.sqrt((a @ a) * (b @ b)))


def neighbour_pairs(shape):
    """Yield the grid points adjacent along each axis, numbered from 1."""
    strides = [int(np.prod(shape[axis + 1 :])) for axis in range(len(shape))]
    for flat in range(int(np.prod(shape))):
        position = np.unravel_index(flat, shape)
        for axis, stride in enumerate(strides):
            if position[axis] + 1 < shape[axis]:
                yield flat + 1, flat + 1 + stride


def read_set(set_dir, point):
    """Return the displacements of one set, with its cell and count checked."""
    ph = phonopy.load(set_dir / "phonopy_disp.yaml", produce_fc=False, log_level=0)
    if not np.allclose(ph.unitcell.cell, point.cell.cell, atol=1e-8):
        print(f"  {set_dir}: not the basis vectors of grid point {point.index + 1}")
    displacements = np.array(ph.dataset["displacements"])
    if len(displacements) != SNAPSHOTS:
        print(f"  {set_dir}: {len(displacements)} snapshots, not {SNAPSHOTS}")
    return displacements


def check_amplitudes(dataset, temperatures):
    """Compare each set with the reference at the temperature it is named after."""
    u = {}  # (grid point numbered from 1, temperature) -> displacements
    print("grid    T        <u^2> / uu")
    for point in dataset.grid_points:
        index = point.index + 1
        reference = reference_u2(point, temperatures)
        for temperature in temperatures:
            set_dir = TRAIN / f"grid-{index:03d}-{temperature}K"
            u[index, temperature] = read_set(set_dir, point)
            ratio = sample_u2(u[index, temperature]) / reference[temperature]
            print(f"{index:4d} {temperature:5d} K   {ratio:10.3f}")
    return u


def check_neighbours(u, temperatures, pairs):
    """Compare the displacement field of neighbouring grid points."""
    print("\n   T    neighbouring grid points, 1 - correlation")
    for temperature in temperatures:
        values = [
            1.0 - correlation(u[a, temperature], u[b, temperature]) for a, b in pairs
        ]
        print(
            f"{temperature:5d} K   median {np.median(values):.1e}   "
            f"max {max(values):.1e}"
        )


def check_merged(dataset, u, temperatures):
    """Check the size and the interleaving of each merged set."""
    print("\nmerged sets")
    expected = len(temperatures) * SNAPSHOTS
    checked = 0
    missing = []
    for point in dataset.grid_points:
        index = point.index + 1
        merged = TRAIN / f"grid-{index:03d}" / "merged.yaml"
        if not merged.exists():
            missing.append(index)
            continue
        ph = phonopy.load(merged, produce_fc=False, log_level=0)
        d = np.array(ph.dataset["displacements"])
        if len(d) != expected:
            print(f"  {merged}: {len(d)} structures, not {expected}")
        for k, temperature in enumerate(temperatures):
            if not np.allclose(d[k :: len(temperatures)], u[index, temperature]):
                print(
                    f"  {merged}: the {temperature} K structures are not at "
                    f"{k}::{len(temperatures)}"
                )
        checked += 1
    print(
        f"  {checked} of {len(dataset.grid_points)} grid points, "
        f"{expected} structures each"
    )
    if missing:
        print(f"  no merged.yaml at grid points {', '.join(map(str, missing))}")


def main():
    """Check the sets of script 5 and the merged sets of script 6."""
    dataset = read_aniso_qha_dataset(DATASET)
    u = check_amplitudes(dataset, TEMPERATURES)
    check_neighbours(u, TEMPERATURES, list(neighbour_pairs(dataset.grid_shape)))
    check_merged(dataset, u, TEMPERATURES)


if __name__ == "__main__":
    main()

The grid points are ordered row-major over dataset.grid_shape, so neighbour_pairs adds the stride of each axis to the flat index of a grid point. A dataset that is not a tensor grid has grid_shape None, and its pairs have to be given to check_neighbours explicitly.

Appendix: spreading the sweep over a cluster#

The 1,025 runs do not talk to each other, so any split of them is allowed. Two are usual.

One job per grid point runs script 7 with the grid point as its argument: 25 jobs of 41 runs. One job per (grid point, temperature) is 1,025 jobs of one run, which suits a queue with a short time limit. The job is then a template with the grid point and the temperature left as placeholders:

Script 10 – one job per grid point and temperature#

#!/bin/bash
#SBATCH --job-name=sscha-GRID-TEMP
#SBATCH --output=sscha-gGRID-tTEMPK.log
grid=train/grid-GRID

phonopy-mlpsscha "$grid/merged.yaml" \
    --mlp "$grid/polymlp.yaml" \
    -t TEMP \
    --snapshots 2000 \
    --iterations 16 \
    --mesh 200 \
    --random-seed 1000 \
    -o "sscha-gGRID-tTEMPK.hdf5"

The #SBATCH lines are Slurm’s, and the resource lines a site asks for belong beside them. Set a name whatever the scheduler. A job read from standard input takes the submitting command as its name, and a queue listing of 1,025 jobs called sbatch says nothing about which run is which.

train/grid-GRID is a relative path, so the job has to start in the directory the sweep was submitted from. Slurm starts it there. Grid Engine needs #$ -cwd in the template, and PBS a cd $PBS_O_WORKDIR before the command.

sed fills the placeholders in, and sbatch reads the job from its standard input when it is not given a file:

for g in {001..025}; do
    for t in {0..400..10}; do
        sed -e "s/GRID/$g/g" -e "s/TEMP/$t/g" job.sh | sbatch
    done
done

The jobs may be submitted in any order and may finish in any order, and a job may be run again. What a rerun reproduces says why, and where that stops.

Appendix: what a rerun reproduces#

What a rerun reproduces says what a seed fixes, what a NumPy upgrade does not, and why an SSCHA run keeps no record of the supercells it drew. Scripts 7 and 10 call phonopy-mlpsscha once per (grid point, temperature), and each of those runs is seeded.

Script 5 passes SEED + int(temperature) to draw_standard_normals, so each of its temperatures is one reproducible draw. The MLPs are fitted to the structures it drew, so another draw there gives another potential. Script 5 writes its draw to normals-*.npz beside the training sets, and sharing one draw of standard normals says how to read one back.

A seeded run draws the same supercells whenever it is made, so the sweep is safe to spread over a cluster and safe to resubmit.

A rerun writing the same file name replaces that file. A rerun writing another name leaves two run files at one (grid point, temperature). Script 8 then stops and prints which grid point and temperature has two.

A run file deleted by mistake stops script 8 as well, and the message lists the first few (grid point, temperature) left without a file. Running those points again is enough, since a rerun draws the supercells the lost run drew.

Anisotropic thermal expansion from free-energy minimization

Contents

Anisotropic thermal expansion from free-energy minimization#

The free energy#

Overview#

0. The equilibrium reference (`phonopy_disp.yaml`)#

1. Build the static grid (internal energy U)#

2. Compute the phonons (phonon grid)#

3. Build the intermediate dataset#

What the builder reads#

How the ordering is checked#

What the file holds#

4. Run the anisotropic QHA#

Supplying \(F_\mathrm{el}\) ready-made#

5. Variant: temperature-dependent force constants from per-grid-point MLPs#

Where this step fits#

One MLP per grid point#

The thermal distribution#

The training displacements#

What the draw leaves at its defaults#

The descriptor and the amount of training data#

Pinning the ridge penalty across the grid#

Validate the MLP#

Computing the free energies with SSCHA#

The options of an SSCHA run and its error#

Gathering the runs into `fph.hdf5`#

Choosing the transient#

Supplying the free energies to the analysis#

Running the analysis#

Appendix: one ladder of descriptors#

Appendix: the gather script#

Appendix: checking the training displacements#

What each check compares#

The script#

Appendix: spreading the sweep over a cluster#

Appendix: what a rerun reproduces#

Anisotropic thermal expansion from free-energy minimization

Contents

Anisotropic thermal expansion from free-energy minimization#

The free energy#

Overview#

0. The equilibrium reference (phonopy_disp.yaml)#

1. Build the static grid (internal energy U)#

2. Compute the phonons (phonon grid)#

3. Build the intermediate dataset#

What the builder reads#

How the ordering is checked#

What the file holds#

4. Run the anisotropic QHA#

Supplying \(F_\mathrm{el}\) ready-made#

5. Variant: temperature-dependent force constants from per-grid-point MLPs#

Where this step fits#

One MLP per grid point#

The thermal distribution#

The training displacements#

What the draw leaves at its defaults#

Sharing one draw of standard normals#

The descriptor and the amount of training data#

Pinning the ridge penalty across the grid#

Validate the MLP#

Computing the free energies with SSCHA#

The options of an SSCHA run and its error#

Gathering the runs into fph.hdf5#

Choosing the transient#

Supplying the free energies to the analysis#

Running the analysis#

Appendix: one ladder of descriptors#

Appendix: the gather script#

Appendix: checking the training displacements#

What each check compares#

The script#

Appendix: spreading the sweep over a cluster#

Appendix: what a rerun reproduces#

0. The equilibrium reference (`phonopy_disp.yaml`)#

Gathering the runs into `fph.hdf5`#