Implementation strategies in phonopy and phono3py

Atsushi Togo; Laurent Chaput; Terumasa Tadano; Isao Tanaka

Source-linked AI summary

Implementation strategies in phonopy and phono3py

Atsushi Togo, Laurent Chaput, Terumasa Tadano, Isao Tanaka

arXiv:2301.05784v1Jan 13, 2023cond-mat.mtrl-sci

TL;DR

Scientific software development and long-term maintenance require implementation knowledge that formal mathematical descriptions may leave implicit. This review explains how computational methods are implemented across phonopy and phono3py, showing applications beyond phonon calculations and providing practical understanding for developers and expert users.

Problem
Scientific software formulas often omit implementation details, while evolving architectures and distributed development make shared knowledge important for maintaining widely used codes.
Method
The review presents the computational methods and implementation strategies used throughout phonopy and phono3py, including supercell, grid, displacement, unfolding, and anharmonic-phonon techniques.
Results
The review describes a collection of methods as implemented in phonopy and phono3py and shows their implementation strategies as a whole.
Takeaways & Limitations
The documented methods are intended to support deeper understanding and development of phonopy and phono3py, while some techniques may also apply to electronic-structure calculations.

Abstract

from arXiv · show

Scientific simulation codes are public property sustained by the community. Modern technology allows anyone to join scientific software projects, from anywhere, remotely via the internet. The phonopy and phono3py codes are widely used open source phonon calculation codes. This review describes a collection of computational methods and techniques as implemented in these codes and shows their implementation strategies as a whole, aiming to be useful for the community. Some of the techniques presented here are not limited to phonon calculations and may therefore be useful in other area of condensed matter physics.

I. INTRODUCTION

This section introduces phonopy and phono3py as open-source phonon-calculation codes and reviews their computational methods and implementation strategies. It emphasizes documentation and community development, then outlines the crystallographic, supercell, force-constant, reciprocal-space, sampling, integration, displacement, and unfolding techniques covered.

Documentation shares software-design choices among users and developers as codes adapt to changing computer architectures, processor cores, memory, and parallel-performance requirements.
Phonopy and phono3py calculate phonon spectra, dynamical structure factors, lattice thermal conductivity, and other physical properties using a supercell-based computational method.
The review provides implementation-level information to help developers understand the codes in depth and participate in their development, while also serving expert users.
The review covers crystal structure and symmetry, reciprocal space, phonon coordinates, supercell geometry, and transformations between force constants and dynamical matrices.
It also describes long-range dipole-dipole corrections, reciprocal-space grids and symmetry, the linear tetrahedron method, finite-temperature random displacements, and phonon band unfolding.

D. Reciprocal space

Reciprocal space provides the framework for describing phonons with wave vectors, reciprocal lattices, symmetry-related wave vectors, and phonon eigenmodes. The section defines the Brillouin zone computationally and addresses its use in phonon calculations and three-phonon scattering.

Phonons are described in reciprocal space using wave vectors, with reciprocal-lattice basis vectors and lattice points constructed from the direct-space primitive cell.
Restricting reciprocal coordinates q_i to [0, 1) defines a primitive cell of reciprocal space.
The star of q is the set {S−⊺q} generated by applying all crystallographic point-group operations to q.
Phonon displacements are expanded in phonon eigenvectors, defining phonon coordinates Q(qν), while frequencies and eigenvectors solve the dynamical-matrix eigenvalue equation.
In phonopy and phono3py, the Brillouin zone is the reciprocal-lattice Wigner–Seitz cell, located by searching reciprocal translations that minimize |q + G|.
For three-phonon scattering, translationally equivalent points on the Brillouin-zone surface are chosen to minimize |q + q′ + q′′| when one or more points lies on that surface.

III. GEOMETRY OF SUPERCELL MODEL

The supercell model is constructed from the primitive cell through an integer matrix, with lattice-point enumeration enabled by matrix reduction and coordinate transformations. Two construction algorithms are described: an auxiliary-supercell method used by default for backward compatibility and a faster alternative for very large supercells.

The supercell basis vectors are represented as the image of primitive-cell basis vectors through an integer matrix Mp→s.
The integer matrix Mp→s is reduced to a diagonal matrix D = P Mp→sQ using unimodular matrices, producing new primitive and supercell vectors with the same volumes and lattices.Because D is diagonal, the new primitive and supercell lattice vectors are collinear, simplifying lattice-point generation.
Lattice points inside the new supercell are generated from D = diag(n1, n2, n3) and mapped into the original supercell using supercell lattice translations and the (mod 1) operation.
An auxiliary supercell uses a diagonal integer matrix Mp→aux = diag(n1, n2, n3), after which its lattice points are shifted into the original supercell and duplicate points are removed.Each original-supercell lattice point can correspond to multiple auxiliary-supercell points; only one is retained.
The auxiliary-supercell construction is the default because it preserves atom indices for backward compatibility, although it is slower than the alternative for very large supercells.

B. Commensurate points

Commensurate reciprocal-space points are those whose phase factors are compatible with supercell translations, enabling force-constant transformations and dynamical-matrix calculations. The implementation handles supercell periodicity by selecting shortest translation vectors and adds a non-analytical correction at general q points to represent long-range dipole-dipole interactions.

Commensurate q points are integer combinations of reciprocal supercell basis vectors satisfying e^iq·R_L = 1 for every supercell lattice vector R_L.For practical calculations, commensurate points within the reciprocal primitive cell are generated using approaches analogous to real-space supercell-point generation.
Supercell force constants are known only within the supercell and are affected there by periodic repetitions of that supercell.
For non-commensurate q, the phase factor depends on which periodic image contains the atom, so the implementation chooses the shortest vectors of R_0κ.Multiple supercell lattice translations may satisfy the atom-pair relation, and the results are stored for repeated use.
The non-analytical term correction adds the reciprocal-space dipole-dipole contribution to the dynamical matrix at general q points, where it is not already included in supercell force constants.For wurtzite-type AlN, the correction is significant near Γ and shows directional dependence toward K and A.
Force constants can be transformed between supercell shapes by oversampling commensurate q points, including embedding anharmonic contributions from a smaller supercell into harmonic force constants of a larger supercell.

VI. REGULAR GRID IN RECIPROCAL PRIMITIVE CELL

The codes discretize reciprocal space with traditional or generalized regular grids, defining microzones and indexing grid points through integer-matrix transformations.

Reciprocal-space integration uses uniformly sampled regular grids, including traditional grids and generalized grids based on conventional-cell reciprocal vectors.
Microzone grid points are generated from integer combinations of microzone basis vectors, optional shifts, and divisions along reciprocal basis directions.
Grid shifts and divisions are chosen so the shifted microzone lattice remains invariant under the crystallographic point group.
A generalized grid uses an integer matrix Mg; det(Mg) gives the number of translationally nonequivalent grid points in the reciprocal primitive cell.
Reciprocal-cell placement may require lattice translations, and periodicity makes distinct coordinate triplets represent equivalent q points.

D. Symmetry of generalized regular grids

Generalized-grid symmetry is tested through integer transformations of grid coordinates and shifts, while BZ-surface points can have multiple equivalent representations.

A generalized grid follows crystallographic symmetry when transformed grid coordinates remain integer under every reciprocal-space point-group operation.
Grid shifts can break symmetry, so the shift must remain invariant modulo 1 under crystallographic point-group operations.
For implementation, shifts normally take values 0 or 1/2 and are doubled so grid arithmetic can use integers instead of floating-point values.
BZ grid points are represented by integer triplets and reciprocal translations chosen to minimize the norm of the represented q point.
Multiple BZ representations can occur when a grid point lies on the BZ surface, so equivalent representations are stored and selected during implementation.

G. Irreducible grid points

The codes exploit crystallographic symmetry to reduce regular-grid calculations, and phono3py further organizes q-point triplets while enforcing momentum conservation.

Irreducible q points represent symmetry stars in the BZ and can reduce computation and memory for single-q phonon properties.
Phono3py identifies irreducible grid points by applying all point-group operations, recovering integer indices, and selecting the minimum valid index.
Dense three-phonon calculations can require 3009 ∼ 10^22 q-point triplets, making direct enumeration impractical.
Three-phonon interactions obey q + q′ + q′′ = G, allowing triplet enumeration to be reduced to pairs of q points on a Γ-centred grid.
Lattice thermal conductivity can be iterated over irreducible q points because heat capacity, lifetime, and group velocity follow phonon-band symmetry.
Triplet sampling fixes q in the irreducible BZ, samples q′ using the subgroup that leaves q invariant, and determines q′′ from momentum conservation.

VII. TETRAHEDRON METHOD

The implemented linear tetrahedron method partitions reciprocal-space microzones into tetrahedra, analytically integrates linearly interpolated functions, and reorganizes results as grid-point weights.

The implementation combines techniques from MacDonald et al. and Blöchl et al. in a modular routine that returns integration weights.
The codes divide each regular-grid microzone into six equal-volume tetrahedra and sum their contributions to obtain the BZ integral.
Linear interpolation within each tetrahedron uses its four vertex values, allowing the tetrahedral integral to be evaluated analytically.
Integration weights are rearranged from tetrahedra to grid points, making the method easier to use and allowing regular-grid symmetry to act directly on the weights.
Weighted densities of states add vertex-function interpolation to tetrahedron integration, reducing to the density of states when F = 1.
Tetrahedron contributions depend on the ordering of the four vertex frequencies, with separate formulae for different frequency intervals and plane-cut geometries.

C. Integration weights

The linear tetrahedron integral is rearranged from tetrahedron contributions into weighted sums over grid points, allowing symmetry handling and shared data structures with smearing methods. For regular grids, local tetrahedron configurations and shifts are predetermined, while symmetry reduction is an approximation when rotated tetrahedra do not map exactly.

Grid-point weights: The integral I(ω) is rearranged from a sum over 6N tetrahedra and four vertices into a weighted sum over grid points.Each grid-point value Fp receives an integration weight wp.
Symmetry and implementation: The rearrangement hides tetrahedron handling, shares data structures with smearing methods, and makes regular-grid symmetry straightforward to apply.Symmetry prepares each grid point’s input before the tetrahedron module processes it.
Grid-point weights: In three dimensions, 24 tetrahedra share each grid point, and the Kronecker selection leaves 24 contributing terms.The 2D illustration has six contributing triangles around a grid point.
Implementation: For regular grids, the four main-diagonal choices and 24 neighboring tetrahedra are hard-coded as a 4×24×4×3 integer array with entries in {−1, 0, 1}.Neighbor positions are represented by shifts from the focused grid point.
Symmetry and implementation: The symmetry relation for integration weights is an approximation because a rotated tetrahedron may not map onto one of the 24 local tetrahedra.The review describes the approximation as very good.
Symmetry and implementation: When the function values obey the grid symmetry, the weighted sum can be computed using irreducible grid points and their multiplicities.The mapping preserves phonon frequencies under the stated symmetry condition.

E. Triplets integration weights

The section extends linear-tetrahedron integration to phonon triplets and explains related implementation strategies for thermal random displacements and defect band unfolding. These procedures retain explicit assumptions about stability, supercell correspondence, and approximate periodic representations.

Triplet integration: For fixed q, momentum conservation selects a single q′′ for each q′, eliminating the q′′ summation in phonon-triplet integration.The resulting quantity has the same functional form as the previously treated integral I(ω).
Triplet integration: Triplet integration weights are computed from function and frequency values at neighboring q′ vertices using the same linear-tetrahedron procedure.The corresponding q′′ neighbor shifts are the negatives of the q′ shifts.
Random displacement generation: Random finite-temperature displacements are generated by sampling real harmonic-oscillator coordinates and using them to construct supercell atomic configurations.The coordinates follow Gaussian distributions associated with Bose–Einstein occupations.
Random displacement generation: Averaging quantities over generated random structures is used to mimic temperature effects on properties such as electronic structure and magnetism.The quantity is computed separately for each structure before averaging.
Random displacement generation: The random-displacement approach is physically valid only when all phonon modes are dynamically stable at every commensurate q point.Harmonic instabilities can occur in some strongly anharmonic systems and high-temperature metastable phases.
Phonon band unfolding: Defect band unfolding computes defective-supercell eigenvectors, maps them through perfect-supercell translations, and decomposes primitive-cell q points as q = ˜q + q⋆.The method assumes one-to-one correspondence between defective-supercell atoms or vacancies and perfect-supercell atoms.
Phonon band unfolding: In the vacancy example, unfolded frequencies are perturbed but the unfolded band structure roughly follows that of the perfect supercell.Only a fraction of modes unfold along the selected L–Γ–L path because the selected q⋆ points are a subset of the commensurate points.

X. CONCLUSION

The review documents phonopy and phono3py implementation details so developers can reconstruct computational methods whose mathematical descriptions may omit practical knowledge. It also identifies regular-grid and linear-tetrahedron techniques as relevant beyond phonon calculations.

Purpose and scope: The review keeps its formalism close to implementation to provide computational-method details useful for distributed scientific software development.The stated audience includes scientific software developers and expert code users.
Purpose and scope: Some reviewed methods are not specific to phonon calculations or the supercell approach.The conclusion explicitly broadens the scope beyond phonon-specific techniques.
Broader applicability: Regular-grid handling and the linear tetrahedron method are identified as applicable to electronic-structure calculations.The review presents these methods as potentially useful to condensed-matter software developers.

Appendix A: Functions used in the linear tetrahedron method

The appendix supplies explicit formulas for the auxiliary functions and weights used by the linear tetrahedron method. It derives these quantities by linearly interpolating within tetrahedra and evaluating occupied subvolumes across frequency-ordering cases.

Definitions: The appendix defines the functions gi, ni, Ii, and Ji used in the linear tetrahedron method.It also records relations among these functions.
Weight formulas: The appendix gives explicit expressions for vertex weights across cases determined by the frequency ordering ωi1 < ωi2 < ωi3 < ωi4.Different occupied-volume geometries produce separate formulas for the weights.
Linear interpolation: Within each tetrahedron, the integrated function is approximated from linearly interpolated vertex values and the tetrahedron’s occupied volume.The contribution uses the tetrahedron volume and the interpolated values at its vertices.
Geometric cases: For partially occupied tetrahedra, the occupied region is decomposed into smaller tetrahedra or represented as the full tetrahedron minus an excluded tetrahedron.The appendix lists the corresponding volumes and weight expressions.
Geometric cases: The listed formulas cover limiting cases, including fully occupied tetrahedra and relations among the auxiliary quantities.The appendix states that a fully occupied tetrahedron is treated directly.

Appendix B: Summation formulae

Appendix B establishes summation formulae for primitive cells and supercells, first in collinear settings and then for arbitrary integer-matrix relations. Noncollinear cases are handled through transformations to collinear cells, while projection periodicity is used in Appendix C.

Appendix B: Summation formulae: When supercell and primitive-cell lattice vectors are collinear, summation formulae (B7) and (B8) are directly valid.This is the starting case for the summation construction.
Appendix B: Summation formulae: For diagonal supercell matrices Mp→s = diag(n1, n2, n3), lattice vectors and conjugated wave vectors fulfill the summation formulae.The formulae allow any reciprocal lattice vector G and supercell lattice vector RL.
Appendix B: Summation formulae: For noncollinear cells, the sums do not easily factor into three geometric series, so validity of the preceding results is not immediately evident.The noncollinear lattice and wave-vector sets are defined within the supercell and reciprocal cell.
Appendix B: Summation formulae: SNF-like transformations can always produce primitive and supercells that are collinear, where the summation results apply.This provides the stated route for treating noncollinear supercell geometries.
Appendix B: Summation formulae: For arbitrary integer matrices Mp→s, unimodular transformations preserve the lattices and yield Eqs. (B7) and (B8).The derivation uses q · RL, G · RL, and G · Rl as multiples of 2π together with det(Mp→s) = det(D).
Appendix B: Summation formulae: In the general construction, wave vectors follow supercell reciprocal vectors within the primitive reciprocal cell, while lattice vectors lie along primitive-cell directions.The sets are defined using integer combinations constrained by the supercell and reciprocal-cell boundaries.
Appendix C: Projection into primitive cell Bloch states: In the projection into primitive-cell Bloch states, a summed position may be rewritten as a position inside the supercell plus a supercell lattice vector.Supercell periodicity of Wlκ(˜q˜ν) and eiq⋆·RL = 1 then support the resulting expression.