(aniso-thermal-expansion)= # 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 . {math}`a`, {math}`b` and {math}`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 ({math}`a`), two for hexagonal, tetragonal and rhombohedral ({math}`a, c`), and three for orthorhombic ({math}`a, b, c`). Cell angles are held fixed, so monoclinic and triclinic crystals are out of scope. This page uses {math}`(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 {math}`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 {math}`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 {math}`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 ```{math} F(a, c; T) = U(a, c) + F_\mathrm{ph}(a, c; T) + F_\mathrm{el}(a, c; T), ``` where the electronic term {math}`F_\mathrm{el}` is optional. {math}`U(a, c)` is the static internal energy, computed by the calculator on the static grid of step 1. {math}`F_\mathrm{ph}(a, c; T)` is the phonon free energy of the harmonic crystal, ```{math} 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 {math}`\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. {math}`a` and {math}`c` are continuous above: {math}`F` is written for any lattice. The calculator returns values only at the sampled cells {math}`(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 {math}`a(T)` and {math}`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: ```{mermaid} flowchart TD EQ(["Equilibrium cell
(phonopy_disp.yaml)"]) EQ --> SC["phonopy-strain-cells
(a, c grid)"] SC --> RELAX{{"calculator relax + static"}} RELAX --> SGRID(["static-grid/grid-NNN
U, F_el"]) SGRID --> PD["generate displacements
per relaxed cell"] PD --> CALCF{{"calculator forces"}} CALCF --> PGRID(["phonon-grid/grid-NNN
disp-*"]) SGRID --> BUILD["phonopy-anisotropic-qha-dataset
--static ... --phonon ..."] PGRID --> BUILD BUILD --> DS(["aniso_qha_dataset.hdf5
cells, U, F_el,
displacements, forces"]) DS --> ANA["phonopy-anisotropic-qha"] ANA --> RES(["a(T), c(T),
alpha_a, alpha_c,
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: ```bash % 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 {math}`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 {ref}`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 -- {math}`a, c` with {math}`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 {math}`U(a_i, c_i)` and, optionally, the electronic states for {math}`F_\mathrm{el}`. ```bash # 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 {math}`c/a`, is the volume path that `phonopy-anisotropic-qha --compare-eos` fits. Equal fractional ranges and equal counts keep {math}`c/a` constant along it, which is the cleanest input to the cross-check. Unequal ranges or counts still give a path, but {math}`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: ```text 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 {math}`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: 1. 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. 2. 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 {math}`F_\mathrm{el}`. A run that writes only `vasprun.xml` carries no eigenvalues, and the dataset is then built without {math}`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 {math}`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. 3. 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. ```{code-block} python :caption: 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 {math}`F_i(T)` of the grid and differentiates it. Either setting biases {math}`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. ```{code-block} python :caption: 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. (anisotropic-qha-build)= ## 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. ``` ```bash % 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: ```bash % 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: ```{code-block} python :caption: 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. (anisotropic-qha-builder-reads)= ### What the builder reads For each grid point the builder reads: - the static single point, giving the internal energy {math}`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 {math}`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 {math}`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 {ref}`phonopy-init -f ` has already collected. This second form is the simpler route when the calculations were not laid out by Script 2: ```bash % 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 ({ref}`--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. (anisotropic-qha-ordering)= ### 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 {math}`U` of one lattice would be combined with the forces of another. - A static single point run on a supercell. Its {math}`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. (anisotropic-qha-file-contents)= ### 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 {math}`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`: ```text / 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: ```bash % 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 {math}`\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 {ref}`exclude_gamma_acoustic_tag`. 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 {math}`U` / {math}`F_\mathrm{ph}` / {math}`F_\mathrm{el}` / total panels, which say which term makes the valley and which one moves it: {math}`U` carries almost the whole curvature, while {math}`F_\mathrm{ph}` and {math}`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 {math}`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 {math}`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 {math}`c/a`, so that `--eos-index` can name a path by hand. Cells that share a {math}`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 {math}`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 {math}`F_\mathrm{el}`. Two meshes can agree closely on {math}`F_\mathrm{el}` and still differ severalfold in {math}`\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: ```{code-block} python :caption: 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: 1. Add the terms at every sampled cell, {math}`F_i(T) = U(a_i, c_i) + F_\mathrm{ph}(a_i, c_i; T)`, with {math}`F_\mathrm{el}(a_i, c_i; T)` when it is included and a {math}`pV_i` term when a pressure is given. 2. Fit a polynomial of total degree {math}`n` in the free lattice parameters to those {math}`F_i(T)`, by least squares. `--polynomial-degree` sets {math}`n`, 3 by default. For {math}`d` free DOF the polynomial has {math}`\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. 3. Minimize that polynomial, from the centroid of the sampled cells and using its analytic gradient. The minimizer is {math}`a(T)`, {math}`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, ```{math} \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 {math}`a(T)`, {math}`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`; {ref}`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 {math}`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 {math}`U` was computed on, and multiplies {math}`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: ```text 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. {math}`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 {math}`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 {math}`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. ```{code-block} python :caption: Script 4 -- {math}`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: ```python 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. ```bash % 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. (anisotropic-qha-temperature-dependent)= ## 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 ({ref}`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 {math}`F(a, c; T)` over the lattice at each temperature. What comes out is {math}`a(T)`, {math}`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, ```{mermaid} flowchart TD DS(["aniso_qha_dataset.hdf5
cells, U, F_el,
harmonic force constants"]) DS --> DISP["thermal displacements
one shared draw of
standard normals"] DISP --> CALCT{{"calculator forces"}} CALCT --> TR(["merged.yaml
per grid point"]) TR --> DEV["train one MLP
per grid point"] DEV --> MLP(["polymlp.yaml
per grid point"]) MLP --> SSCHA["SSCHA
per grid point and temperature"] TR --> SSCHA SSCHA --> FE(["F_ph(T) per
grid point"]) ``` and the analysis they meet in: ```{mermaid} flowchart TD DS2(["aniso_qha_dataset.hdf5"]) FE2(["F_ph(T) per
grid point"]) DS2 --> AQ["run_anisotropic_qha"] FE2 -->|"phonon_free_energies"| AQ AQ --> RES(["a(T), c(T),
alpha_a, alpha_c"]) ``` ### One MLP per grid point The MLP gives the temperature-dependent force constants. {ref}`The SSCHA step ` below computes the anharmonic phonon free energy {math}`F_\mathrm{ph}(a, c; T)` from them. {math}`U(a, c)` comes from the calculator on the static grid of step 1, not from the MLP. The analysis minimizes {math}`F` over the lattice at each temperature. The shape of {math}`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 {math}`U` moves the minimum, and the axial expansions follow that error. That is why {math}`U` is taken from the calculator rather than from the MLP. The same sensitivity is behind the choice of the linear tetrahedron method for {math}`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 {math}`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 {ref}`anisotropic-qha-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. {ref}`The thermal distribution ` writes down the variance {math}`\sigma_{\mathbf{q}\nu}^2` of each mode and the displacements that follow from one draw of standard normals {math}`\xi_{\mathbf{q}\nu}`. A grid point's own frequencies and eigenvectors enter there, which is what {ref}`anisotropic-qha-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. 1. 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. 2. 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 {ref}`anisotropic-qha-normals`. 3. Run the calculator on every supercell, and collect the forces of each set with {ref}`phonopy-init -f `. This writes one `phonopy_params.yaml` per set, holding its displacements, forces and supercell energies. 4. 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. ```{code-block} python :caption: 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: ```bash % 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. {ref}`Merging the temperatures ` gives the reason, which is how `ntrain` and `ntest` cut the merged list: ```{code-block} python :caption: 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 {ref}`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: ```bash % 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 {ref}`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 {math}`|\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. (anisotropic-qha-normals)= ### Sharing one draw of standard normals Drawing the {math}`\xi` once and using the same values at every grid point gives one realization of the randomness over the whole grid. The surface fit then carries it as a smooth function of the lattice parameters. The analysis differentiates that surface, so smoothness matters. `draw_standard_normals` returns the drawn {math}`\xi`, not displacements, and `run(standard_normals=...)` takes a set back. The script above does exactly this. The same {math}`\xi` still gives different displacements at each grid point. Each grid point scales them by its own frequencies and eigenvectors, which makes the draw thermal there. The displacement fields of neighboring grid points then resemble one another. Script 5 writes one `normals-*.npz` per temperature beside the training sets. {ref}`Reproducing and extending the draw ` says what a seed fixes, what a NumPy upgrade does not, and how a saved draw is read back or extended. (anisotropic-qha-validate)= ### The descriptor and the amount of training data {ref}`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. {ref}`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: ```bash % 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 {ref}`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. (anisotropic-qha-sscha)= ### Computing the free energies with SSCHA {math}`F_\mathrm{ph}` is the SSCHA free energy defined in {ref}`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 {ref}`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. {math}`a(T)` and {math}`c(T)` are obtained by minimizing the free-energy surface at each temperature. The axial expansions then come from an Einstein fit of {math}`a(T)` and {math}`c(T)` over the whole temperature range, described in {ref}`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: ```{code-block} bash :caption: 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 {ref}`choosing the transient ` reads off the listing. {ref}`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 {ref}`The options of a run ` covers what each option sets, and {ref}`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. (anisotropic-qha-gather)= ### 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. {ref}`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; {ref}`what a rerun reproduces ` says where that stops. ```bash % ./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`: ```bash % 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. (anisotropic-qha-transient)= ### 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. {ref}`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: ```bash % 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 {math}`F_\mathrm{ph}` from them: ```python run_anisotropic_qha(phonopys, temperatures, internal_energies=internal_energies) ``` Here they are computed elsewhere and passed to `run_anisotropic_qha` instead: ```python 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: ```bash % 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, {math}`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. (anisotropic-qha-sscha-analysis)= ### Running the analysis The file is given to the analysis with `--phonon-free-energies`: ```bash % 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 {math}`da/dT`, {math}`db/dT` and {math}`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 ```{math} a(T) = a_0 + \sum_i A_i \frac{\theta_i}{e^{\theta_i / T} - 1}, ``` with amplitudes {math}`A_i` of opposite sign and Einstein temperatures {math}`\theta_i`, and the same for {math}`b(T)` and {math}`c(T)`. Each term is zero at {math}`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 ```{math} \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`: ```python 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 {math}`a`, {math}`b`, {math}`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 {math}`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: ```{code-block} python :caption: 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. (anisotropic-qha-descriptor-example)= ## Appendix: one ladder of descriptors {ref}`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. (anisotropic-qha-gather-script)= ## Appendix: the gather script What {ref}`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. ```{code-block} python :caption: 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() ``` (anisotropic-qha-check-displacements)= ## 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. {ref}`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 {math}`\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 {math}`g` and {math}`g'` draw from {math}`\tilde{\rho}_{\Phi^{(g)}}(T)` and {math}`\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 ({ref}`anisotropic-qha-normals`), so at a fixed temperature and snapshot their displacements differ only through {math}`\Phi^{(g)}` and {math}`\Phi^{(g')}`. `correlation` measures how close the two sets are as ```{math} \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 {math}`u_{l\kappa j}^{(g,n)}` is snapshot {math}`n` of the set at grid point {math}`g`, and every sum runs over the {math}`N` snapshots and the supercell. Shared normals put {math}`\rho` near 1, and the spacing of the grid sets how far below 1 it falls. `check_neighbours` prints {math}`1 - \rho`. Two draws that do not share their normals are uncorrelated, so {math}`\rho` is 0. A value of {math}`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 ```{code-block} python :caption: 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 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 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 / 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. (anisotropic-qha-distributed)= ## 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: ```{code-block} bash :caption: 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: ```bash 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. {ref}`What a rerun reproduces ` says why, and where that stops. (anisotropic-qha-reproducing)= ## Appendix: what a rerun reproduces {ref}`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 {ref}`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.