# Fileset

[Togo_2023_J._Phys.__Condens._Matter_35_353001.pdf](https://mdr.nims.go.jp/filesets/dfacbe10-eac6-4dd1-8fc6-a8133965f483/download)

## Creator

[Atsushi Togo](https://orcid.org/0000-0001-8393-9766), Laurent Chaput, [Terumasa Tadano](https://orcid.org/0000-0002-8132-2161), Isao Tanaka

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Implementation strategies in phonopy and phono3py](https://mdr.nims.go.jp/datasets/16075f01-2841-4e60-a32b-abf75327c853)

## Fulltext

Implementation strategies in phonopy and phono3pyJournal of Physics: Condensed MatterTOPICAL REVIEW • OPEN ACCESSImplementation strategies in phonopy andphono3pyTo cite this article: Atsushi Togo et al 2023 J. Phys.: Condens. Matter 35 353001 View the article online for updates and enhancements.You may also likeSimulating lattice thermal conductivity insemiconducting materials using high-dimensional neural network potentialEmi Minamitani, Masayoshi Ogura andSatoshi Watanabe-Thermoelectric transport properties ofXAgP (X = Sr and Ba) from first principlesRakshanda Dhawan, Mohd Zeeshan,Tashi Nautiyal et al.-Implications of phonon anisotropy onthermal conductivity of fluorite oxidesSaqeeb Adnan, Miaomiao Jin, Matthew SBryan et al.-This content was downloaded from IP address 144.213.253.16 on 07/08/2023 at 04:19https://doi.org/10.1088/1361-648X/acd831https://iopscience.iop.org/article/10.7567/1882-0786/ab36bchttps://iopscience.iop.org/article/10.7567/1882-0786/ab36bchttps://iopscience.iop.org/article/10.7567/1882-0786/ab36bchttps://iopscience.iop.org/article/10.1088/1361-648X/ac8f7fhttps://iopscience.iop.org/article/10.1088/1361-648X/ac8f7fhttps://iopscience.iop.org/article/10.1088/1361-648X/acd5a1https://iopscience.iop.org/article/10.1088/1361-648X/acd5a1Journal of Physics: Condensed MatterJ. Phys.: Condens. Matter 35 (2023) 353001 (22pp) https://doi.org/10.1088/1361-648X/acd831Topical ReviewImplementation strategies in phonopyand phono3pyAtsushi Togo1,2,∗, Laurent Chaput3, Terumasa Tadano4 and Isao Tanaka2,5,61 Research and Services Division of Materials Data and Integrated System, National Institute forMaterials Science, Tsukuba, Ibaraki 305-0047, Japan2 Center for Elements Strategy Initiative for Structural Materials, Kyoto University, Sakyo, Kyoto606-8501, Japan3 Université de Lorraine, LEMTA, Centre National de la Recherche Scientifique, Unité Mixte deRecherche 7563, BP 70239, 54506 Vandœuvre Cedex, France4 Research Center for Magnetic and Spintronic Materials, National Institute for Materials Science,Tsukuba, Ibaraki 305-0047, Japan5 Department of Materials Science and Engineering, Kyoto University, Sakyo, Kyoto 606-8501, Japan6 Nanostructures Research Laboratory, Japan Fine Ceramics Center, Atsuta, Nagoya 456-8587, JapanE-mail: togo.atsushi@gmail.comReceived 13 January 2023, revised 17 April 2023Accepted for publication 23 May 2023Published 2 June 2023AbstractScientific simulation codes are public property sustained by the community. Modern technologyallows anyone to join scientific software projects, from anywhere, remotely via the internet. Thephonopy and phono3py codes are widely used open-source phonon calculation codes. Thisreview describes a collection of computational methods and techniques implemented in thesecodes and shows their implementation strategies as a whole, aiming to be useful for thecommunity. Some of the techniques presented here are not limited to phonon calculations andmay therefore be useful in other areas of condensed matter physics.Keywords: phonon calculation, computer simulation, software development(Some figures may appear in colour only in the online journal)1. IntroductionSupported by the exponential growth of computer power,computer simulations and scientific software development areessential for condensed matter physics and materials science.Results obtained from computer simulations are getting morerealistic and scientific software development more robust.∗Author to whom any correspondence should be addressed.Original Content from this work may be used under theterms of the Creative Commons Attribution 4.0 licence. Anyfurther distribution of this work must maintain attribution to the author(s) andthe title of the work, journal citation and DOI.Since many scientific codes are distributed under open-source licenses, it is easy to start using scientific softwareand perform computations. Participation in scientific softwaredevelopment has also become easier with the availability ofopen-source compilers and source code editors. Collaborativedevelopment by members at different worksites or locations,i.e. distributed development, has become quite popular for sci-entific software these days. Anyone can contribute to open-source software projects, and the contributions are usuallyunited on internet hosting services for software development.Documentation of scientific software is an important part ofsoftware development, not only for users but also for scientificsoftware developers. For users, it explains how to use the soft-ware. For developers, it describes how to read and write thecode.1361-648X/23/353001+22$33.00 Printed in the UK 1 © 2023 The Author(s). Published by IOP Publishing Ltdhttps://doi.org/10.1088/1361-648X/acd831https://orcid.org/0000-0001-8393-9766https://orcid.org/0000-0002-8132-2161mailto:togo.atsushi@gmail.comhttp://crossmark.crossref.org/dialog/?doi=10.1088/1361-648X/acd831&domain=pdf&date_stamp=2023-6-2https://creativecommons.org/licenses/by/4.0/J. Phys.: Condens. Matter 35 (2023) 353001 Topical ReviewDevelopment of key software can be a project that lasts formany years. Once a key code is developed, it may be used foran unexpectedly long time. However, computer architecturesare constantly evolving. The way to achieve high-performancecomputing may therefore change drastically in a relativelyshort time, on the order of a decade. For example, it is alwaysrequired to follow the increase in the number of cores in pro-cessors. Indeed, algorithms and data structures must be optim-ized for a concurrent use of many cores in a multi-core pro-cessor. Increasing memory space usually allows more flexiblesoftware design. To adapt the code to those changes, the doc-umentation is important to share the meaning of the softwaredesign choices among users and developers.The phonopy and phono3py codes are scientific softwaredeveloped to perform phonon calculations. A variety of phys-ical properties can be calculated using them: the phonon spec-trum, the dynamical structure factor, and the lattice thermalconductivity (LTC), to mention just a few. The computationalmethod is based on the supercell approach. In the computerimplementations, various numerical methods and techniquesare employed. Some of those implemented in the phonopy andphono3py codes are covered in this review. Our motivationfor writing this review is to provide essential information tounderstand the codes in depth and to invite scientific softwaredevelopers in the phonopy and phono3py projects.Although this review is written targeting scientific softwaredevelopers, it may be also useful for expert users of the phon-opy and phono3py codes. The computational methods andtechniques underlying the phonon calculations are describedas implemented to the phonopy and phono3py codes.In section 2, the representation of the crystal structureand the crystal symmetry are presented in a crystallographicway. The reciprocal space of the crystal and the crystal sym-metry in the reciprocal space are also described. Then, thephonon coordinates are introduced. The concept of Brillouinzone (BZ) required for the phonon calculation is brieflydefined.In section 3, the geometry of the supercell structure modelis described. The supercell geometry is associated with theprimitive cell by an integer matrix and how to construct thesupercell geometry using the integer matrix is explained.In section 4, the transformation between the supercell forceconstants and the dynamical matrices is presented. Since thephonopy and phono3py codes employ the supercell approach,the force constants elements are limited to within the super-cell. The transformation, therefore, needs to be performedwithspecial care. Such details are given in this section.In section 5, the treatment of long range dipole–dipoleinteraction in the dynamical matrix, as implemented in thephonopy code, is explained. The implementation is mainlybased on [1], and its application to the supercell approach isreviewed.In section 6, the traditional and generalized regular gridsused to sample reciprocal space discretely are defined. Thesymmetry treatment of those grid points, used in the phonopyand phono3py codes, is also presented.In section 7, a linear tetrahedron method is described usingthe formulation as implemented in the phonopy and phono3pycodes. The implementation is mainly based on [2, 3]. Themethod is then applied to the three-phonon scattering.In section 8, a scheme to generate random atomic displace-ments in the supercell at finite temperatures is presented. Thisis achieved by superpositions of displacements correspondingto randomly sampled harmonic oscillators.In section 9, a band unfolding technique implemented inthe phonopy code following [4] is explained. The formulationspecific to the phonon calculation with the supercell approachis provided.2. Crystal and symmetry2.1. Crystal structureAs shown in figure 1, a crystal model is defined by a col-lection of unit cells periodically repeated in the three direc-tions of space. Each unit cell contains atoms. We normallychoose a conventional set of basis vectors to span the unit cell,(a,b,c), so that it naturally represents the crystal symmetryalong with the choice of the origin of atomic positions [5].Positions of atoms are represented either in Cartesian coordin-ates or in crystallographic coordinates. Denoting the crystallo-graphic coordinates of an atom as x= (x1,x2,x3)⊺, its positionR is given byR= x1a+ x2b+ x3c= (a,b,c)x1x2x3= (a,b,c)x. (1)In the above equation, if the basis vectors are represented usingcolumn vectors containing their Cartesian components, then(a,b,c) becomes a matrix, and R a column vector containingthe Cartesian coordinates of the atomic position vector,RxRyRz=ax bx cxay by cyaz bz czx1x2x3 . (2)In the following we will use both representations. Dependingon the context, (a,b,c) can therefore be considered as a rowvector containing three elements of a vector space, or a 3 by 3matrix.It is assumed that elements of x are in the interval [0,1) todescribe the structure. However it is often the case that we haveto compare two positions of atoms which may be in differentunit cells, possibly after a space group operation. For example,to identify two coordinates x and x ′ which correspond to thesame location within the unit cell, but are possibly shifted bya lattice vector, it is convenient to bring each element of thedifference ∆x= x ′ − x into the interval [−0.5,0.5) and thenconfirm that the length of∆x is smaller than a tolerance ϵ. This2J. Phys.: Condens. Matter 35 (2023) 353001 Topical ReviewFigure 1. Crystal structure model of β-Si3N4. See also figure 2 forthe space group symmetry. (a) Unit cell containing atoms. (b) Basisvectors of the unit cell. (c) Crystal lattice. The filled circle symbolsshow the lattice points. In this model, the unit cell is assumed to beprimitive.is implemented in the phonopy and phono3py codes every-where by rounding components of ∆x to the nearest integer(nint(∆x)) and checking |(a,b,c)[∆x− nint(∆x)]|< ϵ.2.2. Space group operationVarious symmetries of crystals and phonons are used in thephonopy and phono3py codes. The most important one is thespace group symmetry. A space group operation is composedof a (proper or improper) rotation S and a translation τ . It isoften written with a composite notation (S, τ) which sendsposition R toR ′ = (S, τ)R≡ SR+ τ. (3)Using crystallographic coordinates this is written as(a,b,c)x ′ = S(a,b,c)x+(a,b,c)w (4)orx ′ = (a,b,c)−1S(a,b,c)x+w≡ Sx+w, (5)where S can be shown to be a unimodular matrix. w con-tains the components of the translation vector in the crystal-lographic basis. If Cartesian coordinates were used instead,the matrix representing the rotation would be orthogonal.Representing the space group by S= {(S, τ)}, the crystal-lographic point group is given by P= {S|(S, τ) ∈ S}. Anexample of space group operation is shown in figure 2. Thecrystal structure overlaps with itself after a rotation of 60◦along the z axis and a + 12c shift. ThereforeFigure 2. Symmetry operations applied on the unit cell of β-Si3N4whose space-group type is P63/m (No. 176). The left and middlefigures illustrate identity and six-fold screw axis operations. Theirnotations follow the international tables for crystallography volumeA [5]. The right figure shows the crystal structure viewed from thetop (z-axis) that depicts the six-fold screw axis passing through theorigin O.S=1 1̄ 01 0 00 0 1 ,w= 001/2 , (6)must be a space group operation. It means that any atom in theunit cell at point x is sent by the space group operation to anew pointx ′ =1 1̄ 01 0 00 0 1x+ 001/2 . (7)The new point x ′ can be located out of the original unit cell,but an atom with the same atomic type must be found at thisnew location.2.3. Primitive cellIn most cases, a unit cell as described in section 2.1 is usedas the input crystal structure model of a phonon calculation.The most plausible choice of the unit cell is either a primit-ive cell or a conventional unit cell. The former represents theminimum unit of periodicity of a crystal lattice. It is there-fore the one used for Fourier expansions and is necessary fora reciprocal space representation of phonon properties. TheBloch theorem evidenced this lattice translational symmetry.The latter follows the crystallographic convention. Althoughit provides intuitive shapes of the unit cells (cubic, hexagonal,etc), it may contain several lattice points. Figure 3 shows theconventional unit cell and a primitive cell for silicon crystal-lized within the face-centred-cubic structure.The conventional unit cell structure is uniquely defined tofollow the crystallographic convention, [5] within the freedomdue to Euclidean normalizer [5]. On the other hand, the choiceof the primitive cell is not unique. This is inconvenient fora systematic handling of crystal structure models in phononcalculations. Therefore, phonopy suggests a predefined trans-formation matrix Pprim from the conventional unit cell to theprimitive cell. Pprim is used as(ap,bp,cp) = (ac,bc,cc)Pprim, (8)3J. Phys.: Condens. Matter 35 (2023) 353001 Topical ReviewTable 1. Choices of transformation matrices Pprim of equation (8)used in the phonopy and phono3py codes. The subscripts X of thematrices PX indicate the centring types: A, B, C for the base centringtypes, I and F for the body and face centring types, respectively, andR for the (obverse) rhombohedral centring type.PA =1 0 00 121̄20 1212, PB = 12 0 1̄20 1 012 0 12, PC = 1212 01̄212 00 0 1,PI = 1̄21212121̄21212121̄2, PF =0 121212 0 121212 0, PR = 231̄31̄313132̄3131313.Figure 3. Conventional unit cell and primitive cell of silicon. Thespace group type is Fd3̄m (No. 227). Since the conventional unit cellcontains four lattice points with eight atoms, the primitive cell ischosen to have two atoms by det(PF) = 1/4.where (ap,bp,cp) are the basis vectors of the primitive celland (ac,bc,cc) those of the conventional unit cell. The choicesof the transformation matrices used in the phonopy andphono3py codes are shown in table 1. For example, the trans-formation matrix used for silicon in figure 3 is Pprim = PF.These matrices are equivalent to those presented in table 2of [6] although those are given for the reciprocal space basisvectors.The crystal lattice is defined as the set of integer linear com-binations of primitive lattice vectors. A lattice vector Rl cantherefore be written asRl = l1ap + l2bp + l3cp = (ap,bp,cp)l (li ∈ Z). (9)In figure 1, if the unit cell is assumed to be primitive, thenfigure 1(c) shows the crystal lattice. The lattice point maythen be chosen to coincide with the origin of each unit cell.Sometimes it is also useful to regard the crystal lattice as gen-erated using the lattice translations of the space group, (I, t),where I and t denote the identity operation and lattice transla-tion, respectively.2.4. Reciprocal spaceThe crystal structure, and the symmetry explained above, aredefined in direct (real) space. In this section, crystallographyis briefly summarized in reciprocal space. This is quite use-ful for our purpose since phonons are described in reciprocalspace using wave vectors. It is convenient to introduce a recip-rocal lattice whose basis vectors (a∗p ,b∗p ,c∗p) are solutions ofthe equation(ap,bp,cp)⊺(a∗p ,b∗p ,c∗p) = 2π1. (10)This equation can be solved to givea∗p = 2πbp × cpVc, b∗p = 2πcp × apVc, c∗p = 2πap ×bpVc, (11)where Vc = ap · (bp × cp) = bp · (cp × ap) = cp · (ap ×bp).The reciprocal lattice points are then obtained from the integerlinear combinations of reciprocal basis vectors,G= G1a∗p +G2b∗p +G3c∗p = (a∗p ,b∗p ,c∗p)G (Gi ∈ Z). (12)Similarly to the atomic positionsR, wave vectors q are definedfrom real linear combinations of the reciprocal lattice vectors,q= q1a∗p + q2b∗p + q3c∗p = (a∗p ,b∗p ,c∗p)q. (13)If the qi are restricted to the interval [0,1), a primitive cell ofreciprocal space is obtained.If we consider a component of a Fourier expansion, fq(R) =eiq·R, it is transformed by an operation of the crystallographicpoint group to the functionSfq(R) = fq(S−1R) = eiq·S−1R = eiSq·SS−1R = eiSq·R= fSq(R). (14)For this reason, the image of a wave vector through an opera-tion of the crystallographic point group is defined to beq ′ = Sq⇐⇒ (a∗p ,b∗p ,c∗p)q′ = S(a∗p ,b∗p ,c∗p)q (15)orq ′ = (a∗p ,b∗p ,c∗p)−1S(a∗p ,b∗p ,c∗p)q (16)= (ap,bp,cp)⊺S(ap,bp,cp)−⊺q (17)= S−⊺q. (18)The set {S−⊺q} with S in the crystallographic point group{S} is called the star of q.2.5. Phonon coordinatesAtoms vibrate in the vicinity of their equilibrium positions incrystals. Instantaneous and equilibrium positions of atoms incrystals are denoted by Rlκ and R0lκ, respectively, where l andκ are used to label the lattice points and atoms in the primitivecell of l, respectively. Displacements of atoms are written asulκ = Rlκ −R0lκ.In this review, we define the dynamical matrix as [7]Dκα,κ ′α ′(q) =1√mκmκ ′∑l ′Φ0κα,l ′κ ′α ′eiq·(R0l ′κ ′−R00κ).(19)α and mκ denote the index of the Cartesian coordinates andthe mass of the atom κ, respectively. Φlκα,l ′κ ′α ′ are the har-monic force constants, defined as the second derivative of theenergy with respect to atomic positions, and evaluated at the4J. Phys.: Condens. Matter 35 (2023) 353001 Topical Reviewequilibrium positions [7]. Phonon frequency ωqν and eigen-vectorWκα(qν) are obtained as the solution of the eigenvalueequation of the dynamical matrix in equation (19), which iswritten as∑κ ′α ′Dκα,κ ′α ′(q)Wκ ′α ′(qν) = ω2qνWκα(qν). (20)ν labels the phonon band index, and the composite index qνis used to consider a phonon mode.Equation (20) can also be written in matrix form asD(q)W(q) =W(q)Ω2(q), (21)or, because D(q) is hermitian,D(q) =W(q)Ω2(q)W†(q), (22)where Ω2(q) is the diagonal matrix whose diagonalelements are ω2qν . Each column of W(q) contains aneigenvector Wκα(qν) corresponding to a different bandindex ν. Elements of each column are ordered as(Wκ1x,Wκ1y,Wκ1z, . . . ,Wκna x,Wκna y,Wκna z), where na is thenumber of atoms in the primitive cell.Because Dκα,κ ′α ′(q) = D∗κα,κ ′α ′(−q), we haveωqν = ω−qν , (23)and we can chooseWκα(qν) =W∗κα(−qν). (24)Moreover eigenvalues and eigenvectors inherit symmetryproperties of the force constants. It can be shown [8] that undera space group operation (S, τ) the eigenvalues and eigen-vectors transform asωSqν = ωqν , (25)Wκ ′(Sqν) = SWκ(qν)e−iSq·τ , (26)where Wκ(qν) represents the phonon eigenvector of atom κin Cartesian coordinates, and κ ′ is the image of atom κ underthe space group operation.An actual displacement can be described from a linearcombination of phonon eigenvectors. This defines the phononcoordinates Q(qν) asulκα =∑qνQ(qν)1√NmκWκα(qν)eiq·R0lκ (27)=∑qνQ(qν)ulκα(qν), (28)where N is the number of lattice points in crystal. ulκα(qν)is defined for the purpose of convenience in the followingsections. The previous equation can be inverted to giveQ(qν) =∑lκαulκα√mκNW∗κα(qν)e−iq·R0lκ . (29)Figure 4. Brillouin zone (BZ) of β-Si3N4. See also figures 1 and 2for the cell shape in direct space. (a) The first BZ. (b) BZs viewedfrom c∗p axis direction. Symbols of the special points and pathsfollow the Bilbao crystallographic server [6].2.6. BZSymmetry property of phonons in reciprocal space is best rep-resented in the BZ [6, 8–12]. In the phonopy and phono3pycodes, the BZ is defined as a Wigner–Seitz cell of the recip-rocal lattice. To check if a q point belongs to the BZ weproceed as follows. First the basis vectors of the Niggli cell[13–16] are determined. The Niggli cell is a cell with theshortest possible reciprocal basis vectors. Then using integerlinear combinations of those Niggli basis vectors, we searchfor the shortest |q+G|. Finally q+G is used as the q point inthe BZ.Three q points are needed for the three-phonon scatteringconsidered in the phono3py code. In case one or more of thosethree points is on the BZ surface, we choose the translationallyequivalent points on the BZ surface which minimize |q+q ′ +q ′ ′|.As an example, the BZ of β-Si3N4 is presented in figure 4.The basal plane has the hexagonal shape and c∗p is longer thana∗p and b∗p because cp is shorter than ap and bp. By definition,a∗p/2, b∗p/2, and c∗p/2 are located on the BZ surface. The highsymmetry points and paths in the BZ have special symbols asshown in figure 4(b) [6].3. Geometry of supercell model3.1. Supercell constructionIn the phonopy and phono3py codes, the supercell approachis employed. The supercell is defined by multiple primitivecells, so that the basis vectors of the supercell (as,bs,cs) canbe represented as the image of the primitive cell basis vectorsthrough an integer matrixMp→s,(as,bs,cs) = (ap,bp,cp)Mp→s. (30)Like the unit cell, the supercells are arranged on the supercelllattice. The supercell lattice points are labeled by L, and loc-ated from the vectors RL.The supercell model is constructed usingMp→s. The latticepoints within the supercell are used to perform the summationfound in equation (19). These lattice points can be elegantlyobtained using the approach reported by Hart and Forcade5J. Phys.: Condens. Matter 35 (2023) 353001 Topical Review[17]. The integer matrix Mp→s is first reduced to a diagonalinteger matrix D= PMp→sQ by the unimodular matrices PandQ, where we choose det(P) = 1 and det(Q) = 1. This mat-rix decomposition can be the Smith normal form (SNF), how-ever, it is unnecessary to strictly follow the definition of theSNF.Using this transformation, equation (30) can be written(ãs, b̃s, c̃s) = (ãp, b̃p, c̃p)D, (31)with(ãs, b̃s, c̃s) = (as,bs,cs)Q, (32)(ãp, b̃p, c̃p) = (ap,bp,cp)P−1. (33)(ãs, b̃s, c̃s) and (ãp, b̃p, c̃p) define a new supercell and a newprimitive cell, respectively. Because det(P) = det(Q) = 1,they have the same volumes as the former ones and gener-ate the same lattices. The benefit of using those new primitiveand supercell comes from the fact that D is diagonal. The newprimitive and supercell lattice vectors are therefore collinear.We can easily find the lattice vectors located within the newsupercell. If we write D= diag(n1,n2,n3), they are given bythe vectorsRm = (ãp, b̃p, c̃p)m1m2m3= (ãp, b̃p, c̃p)m, (34)with mi ∈ {0,1, . . . ,ni − 1}. This can also be writtenRm = (ãs, b̃s, c̃s)D−1m (35)= (as,bs,cs)QD−1m. (36)Therefore, the lattice points within the original supercell havecoordinatesxs = QD−1m (mod 1) (37)along the supercell lattice vectors (as,bs,cs). The (mod 1)operation is used to shift the lattice vectors from within thenew supercell to within the original supercell.In figure 5, an example of a supercell construction for a two-dimensional lattice is presented.D= PMp→sQ is computed as(2 00 4)=(0 1−1 1)(2 2−2 2)(0 −11 −1). (38)This gives(ãp, b̃p) = (ap,bp)P−1 = (ap +bp,−ap), (39)(ãs, b̃s) = (as,bs)Q= (bs,−as −bs). (40)As can be seen in figure 5, (ãp, b̃p) and (ãs, b̃s) are mutuallyparallel, and the new supercell is simply built by 2× 4 unitsFigure 5. Lattices and basis vectors of primitive cell (ap,bp),supercell (as,bs), new primitive cell (ãp, b̃p), new supercell (ãs, b̃s),and auxiliary supercell (aaux,baux). Circle and small filled circlesymbols depict lattice points in the supercell and in the newsupercell, respectively. The latter points can be brought by thesupercell lattice translations to the former points uniquely.of the new primitive cell. Since the lattices generated by theoriginal or new basis vectors are the same, the lattice pointsinside the new supercell are simply brought into the originalsupercell by supercell lattice translations.Another way of constructing the supercell is to employ anauxiliary supercell,(aaux,baux,caux) = (ap,bp,cp)Mp→aux, (41)which is defined by a diagonal integer matrix Mp→aux =diag(n1,n2,n3). This gives(aaux,baux,caux) = (as,bs,cs)M−1p→sMp→aux. (42)We choose M−1p→sMp→aux to be an integer matrix. To satisfyit, in the implementation, the diagonal elements of Mp→auxare determined by making the smallest parallelepiped of(aaux,baux,caux) that includes the parallelepiped of (as,bs,cs).The lattice points within this auxiliary supercell are thenRm = (ap,bp,cp)m1m2m3= (ap,bp,cp)m, (43)with mi ∈ {0,1, . . . ,ni − 1}. Their coordinates along the basisvectors of the auxiliary supercell arexaux =M−1p→auxm (44)while their coordinates along the basis vectors of the originalsupercell arexs =M−1p→sMp→auxxaux (mod 1). (45)As before, we used the (mod 1) operation to shift the lat-tice points from within the auxiliary supercell to within theoriginal supercell. Notice that every lattice point within theoriginal supercell is obtained from det(M−1p→sMp→aux) lattice6J. Phys.: Condens. Matter 35 (2023) 353001 Topical Reviewpoints within the auxiliary supercell. Only one of them shouldbe conserved.An example is shown in figure 5. Lattice points in the aux-iliary supercell are brought into the supercell by the supercelllattice translations. But since det(M−1p→sMp→aux) = 2, each lat-tice point within the original supercell is obtained from twolattice points in the auxiliary supercell.For backward compatibility, to preserve indices of atomsin the supercell, it is the latter way of the supercell construc-tion which is used as default in the phonopy and phono3pycodes. However, the latter algorithm is much slower than theformer, and therefore the former should be used for very largersupercells.3.2. Commensurate pointsFrom equations (30) and (10), we have(a∗p ,b∗p ,c∗p) = (a∗s ,b∗s ,c∗s )M⊺p→s. (46)(a∗s ,b∗s ,c∗s ) are the reciprocal basis vectors associated with thesupercell. q points given by integer linear combinations ofthe reciprocal supercell basis vectors are called commensuratewith the supercell because they fulfil eiq·RL = 1 for any super-cell lattice vector RL.For practical purposes, we are interested in the commensur-ate points located within the reciprocal primitive cell. Sinceequation (46) has the same form as equation (30), the sameapproaches as those explained in section 3.1 can be used forthe generation of the commensurate points within the recip-rocal primitive cell.4. Transformation between force constants anddynamical matricesIn equations (19) and (20), force constants are transformedto phonon eigenvalues and eigenvectors. In principle, inequation (19) the summation over l′ is performed for all lat-tice vectors. However, in the supercell approach, it is notΦ0κα,l ′κ ′α ′ which is obtained from numerical simulations,but ΦSC0κα,l ′κ ′α ′ =∑L ′Φ0κα,L ′+l ′κ ′α ′ , [18] with l′ restrictedto the supercell. It means that the force constants are onlyknown within a restricted region of the crystal, the supercell,and within this region, they are contaminated by the periodicrepetitions of the supercell. This also means that the supercellforce constants ΦSC0κα,l ′κ ′α ′ are periodic under a supercell lat-tice translation while the phase factor eiq·(R0l ′κ ′−R00κ) is not,except when q is a wave vector commensurate with the super-cell. There is therefore an ambiguity about where to choosethe atom (l ′κ ′). For commensurate wave vectors, it could bechosen in any periodic repetition of the supercell containingatom (0κ) without changing the value of the phase factor.However, for other values of q, this value would be changed,and therefore this choice matters.In the phonopy and phono3py codes, following [18], thephase factors given by the shortest vectors of R0l ′κ ′ +RL−R00κ are chosen, where RL are the positions of the supercelllattice points L. This is implemented by searching L by{L}κκ ′l ′ = argminL(|R0l ′κ ′ +RL−R00κ|), (47)and the results are stored and used many times in the cal-culation. Since multiple L can be found for each pair ofatoms (0κ) in the primitive cell and (l ′κ ′) in the supercellcell, equation (47) gives a set of L. The dynamical matrix ofthe supercell is computed from the supercell force constantsΦSC0κα,l ′κ ′α ′ averaging the phase factors of {L} as follows, [18]DSCκα,κ ′α ′(q) =1√mκmκ ′∑l ′ΦSC0κα,l ′κ ′α ′× 1|{L}κκ ′l ′ |∑{L}κκ ′ l ′eiq·(R0l ′κ ′+RL−R00κ).(48)Equation (19) can be inverted to compute force constantsfrom N dynamical matrices, [19]Φ0κα,l ′κ ′α ′ =√mκmκ ′N∑qDκα,κ ′α ′(q)e−iq·(R0l ′κ ′−R00κ).(49)Dκα,κ ′α ′(q) can be obtained from known eigenvectorsand eigenvalues, as in equation (22), or calculated byequation (48). ΦSC0κα,l ′κ ′α ′ can be obtained from DSCκα,κ ′α ′(q)by restricting the above q sum to the wavevectors q⋆ com-mensurate with the supercell. Equation (49) can also be usedto transform supercell force constants in a different super-cell shape by oversampling q⋆ points. Indeed, we obtainΦSCL←S0κα,l ′Lκ′α ′ =1N∑q⋆Le−iq⋆L ·(R0l ′Lκ′−R00κ)∑l ′SΦSCS0κα,l ′Sκ′α ′1|{LS}κκ ′l ′S |∑{LS}κκ ′ l ′Seiq⋆L ·(R0l ′Sκ′+RLS−R00κ) , (50)where subscripts of variables XL and XS mean those definedin different supercell shapes, e.g. larger (L) and smaller (S)supercells, respectively. The original supercell force constantsΦSCS0κα,l ′Sκ′α ′ are given in the smaller supercell. Commensuratepoints q⋆L are sampled with respect to the larger supercell. Thetransformed supercell force constants ΦSCL←S0κα,l ′Lκ′α ′ are given in7J. Phys.: Condens. Matter 35 (2023) 353001 Topical Reviewthe larger supercell. A possible application is to embed anhar-monic contribution obtained through self-consistent harmonicapproximation using a smaller supercell into the harmonicforce constants of a larger supercell [20].5. Non-analytical term correctionLong-range dipole–dipole interactions are difficult to capturein a supercell approach. Therefore it is treated with the help ofa model, the so-called non-analytical term correction (NAC)[1, 21–23].At the commensurate points, this contribution is alreadyincluded in equation (48) via the supercell force constants.However, at general q points it is not. Gonze and Lee [1]have formulated this contribution to the dynamical matrix,and in the phonopy and phono3py codes only the reciprocalspace term of the dipole–dipole interaction contribution to thedynamical matrix is calculated. It is given byDDDκα,κ ′α ′(q) =4πVce24π ϵ01√mκmκ ′∑ββ ′Z∗κ,βαZ∗κ ′,β ′α ′×∑Q=G+qQβQβ ′∑γγ ′Qγϵγγ ′Qγ ′× eiG·(R00κ−R00κ ′ ) exp−∑γγ ′Qγϵγγ ′Qγ ′4Λ2 ,(51)where β and γ indicate the Cartesian coordinates. If the polar-ization is called P, then Z∗κ,βα = Vce∂Pβ∂R0κα|E=0 are the Borneffective charges at zero electric field, E= 0. ϵγγ ′ is the high-frequency static dielectric constant tensor. Λ is a parameteradjusted together with the cutoff radius used for the summationover G. A translational invariance condition can be applied toequation (51) by [1]DDDκα,κ ′α ′(q)←DDDκα,κ ′α ′(q)− δκκ ′∑κ ′ ′√mκ ′ ′mκDDDκα,κ ′ ′α ′(q= 0).(52)As said already, at the commensurate points, this dipole–dipole contribution is already included in equation (48) via thesupercell force constants. Therefore, equation (51) is used forthe interpolation of the dynamical matrix at general q points.The procedure is reported in [1], which is described shortlyas follows. At the commensurate points q⋆, the short-rangedynamical matrix is calculated asDSR,SCκα,κ ′α ′(q⋆) = DSCκα,κ ′α ′(q⋆)−DDDκα,κ ′α ′(q⋆). (53)Next, the short-range supercell force constantsΦSR,SC0κα,l ′κ ′α ′ areobtained from {DSR,SCκα,κ ′α ′(q⋆)} using equation (49), then theshort-range dynamical matrix DSR,SCκα,κ ′α ′(q) is calculated at thegeneral q point by equation (48), i.e.Figure 6. Calculated phonon band structure of wurtzite-type AlNusing the 5× 5× 3 supercell with (solid curves) and without(dashed curves) non-analytical term correction (NAC). Circle andsquare symbols depict phonon frequencies at q→ 0, where filledand open symbols indicate phonon modes calculated with andwithout NAC, respectively. Wave vector path was selected bySeek-path [24, 25].{DSR,SCκα,κ ′α ′(q⋆)}−−−−→Eq. (49)ΦSR,SC0κα,l ′κ ′α ′ −−−−→Eq. (48)DSR,SCκα,κ ′α ′(q).(54)Finally, the dynamical matrix at the general q point with NACis obtained byDNACκα,κ ′α ′(q) = DSR,SCκα,κ ′α ′(q)+DDDκα,κ ′α ′(q). (55)An example of the application of NAC towurtzite-type AlNis shown in figure 6. As can be seen, the correction is signific-ant near the Γ point. It also shows the directional dependencenear theΓ, in the direction of K andA, due to the non-sphericalsymmetry of ϵγγ ′ and Z∗κ,βα of AlN.6. Regular grid in reciprocal primitive cellWhen phonon properties at q points need to be integratedover the BZ, a technique often used is to discretize the recip-rocal space. Usually, reciprocal primitive cells are uniformlysampled by a regular grid such as [26], even if other choicesare possible [27]. Traditionally, the regular grid is definedby evenly dividing the reciprocal basis vectors, as shownin figure 7(a). However, other regular grids may be chosen,as shown in figure 7(b). For example, the grid shown infigure 7(b) is defined by dividing evenly the reciprocal basisvectors of the conventional unit cell instead. This is a type ofgeneralized regular grid, [28–31] as will be explained later.6.1. Traditional regular gridTraditionally, the volume of the reciprocal primitive cell isdivided into uniform microzones so that the basis vectors of8J. Phys.: Condens. Matter 35 (2023) 353001 Topical ReviewFigure 7. Different Γ-centred regular grids in plane reciprocalprimitive cell (c2mm). Filled circle symbols depict the grid points.Open circle symbols show the grid points equivalent to some othergrid points by periodicity. Shaded area indicates a microzone [2] anda∗m and b∗m give the basis vectors of the microzone. (a) Traditionalregular grid. (b) A type of generalized regular grid [28–31].the reciprocal primitive cell are simply integer multiples ofthe basis vectors of each microzone (a∗m,b∗m,c∗m). The equationwhich defines the microzone basis vectors is then(a∗p ,b∗p ,c∗p) = (a∗m,b∗m,c∗m)D (56)with D= diag(n1,n2,n3), ni ∈ N.A 2D example of this microzone is shown in figure 7(a).The q points of the grid points are represented by integer lin-ear combinations of the microzone basis vectors plus a pos-sible rigid shift (s1,s2,s3)⊺. Therefore their coordinates in the(a∗p ,b∗p ,c∗p) basis areq1q2q3=(m1 + s1)/n1(m2 + s2)/n2(m3 + s3)/n3 , (57)mi ∈ {0,1, . . . ,ni − 1}, 0⩽ si < 1.To conserve symmetry, ni and si are chosen so that the micro-zone lattice with the shift is invariant under the crystallo-graphic point group P.6.2. Generalized regular gridA generalized regular grid is defined using a conventional unitcell related to the primitive cell by equation (8). The reciprocalconventional unit cell is therefore(a∗c ,b∗c ,c∗c ) = (a∗p ,b∗p ,c∗p)P⊺prim. (58)Microzones can be defined considering that their basis vectorsare integer divisions of the reciprocal basis vectors of the con-ventional unit cell. This is written as(a∗c ,b∗c ,c∗c ) = (a∗m,b∗m,c∗m)diag(N1,N2,N3), Ni ∈ N. (59)A 2D example of this microzone is shown in figure 7(b). Fromequations (58) and (59), we have(a∗p ,b∗p ,c∗p) = (a∗m,b∗m,c∗m)diag(N1,N2,N3)P−⊺prim= (a∗m,b∗m,c∗m)Mg. (60)The grid matrixMg = diag(N1,N2,N3)P−⊺prim is an integer mat-rix. det(Mg) is the number of translationally nonequivalentgrid points in the reciprocal primitive cell. When Mg is nota diagonal matrix, the grid generated by equation (60) is ageneralized regular grid, otherwise we obtain the traditionalregular grid as described in section 6.1. Notice however thatthe generalized regular grid may be transformed into the tradi-tional regular grid of some reciprocal cell using an SNF kindof decomposition, as explained in section 3.1. This is shownin the next section.6.3. Indexing of grid pointsThe integer matrix Mg can be reduced to a diagonalinteger matrix such as D= PMgQ as has been employed insection 3.1. This property of the matrix decomposition issimilarly used to index grid points by integer numbers of{0,1, . . . ,det(D)− 1} [17] in the phono3py code. Note thatdet(D) = det(Mg). Equation (60) is rewritten as(a∗p ,b∗p ,c∗p)Q= (a∗m,b∗m,c∗m)P−1D. (61)Denoting(ã∗p , b̃∗p , c̃∗p) = (a∗p ,b∗p ,c∗p)Q, (62)(ã∗m, b̃∗m, c̃∗m) = (a∗m,b∗m,c∗m)P−1, (63)we have(ã∗p , b̃∗p , c̃∗p) = (ã∗m, b̃∗m, c̃∗m)D. (64)Since Q is a unimodular matrix, (ã∗p , b̃∗p , c̃∗p) and (a∗p ,b∗p ,c∗p)generate the same reciprocal primitive lattice. Similarly(ã∗m, b̃∗m, c̃∗m) and (a∗m,b∗m,c∗m) generate the samemicrozone lat-tice due to the unimodular matrix P−1.Equation (64) has the same form as equation (56) withD= diag(n1,n2,n3). Therefore the q points of the grid pointsare calculated similarly as equation (57) but in the basis(ã∗p , b̃∗p , c̃∗p),q= (ã∗p , b̃∗p , c̃∗p)(m1 + s̃1)/n1(m2 + s̃2)/n2(m3 + s̃3)/n3 . (65)Notice that in the above equation we use s̃i rather than si. Infact for practical purposes, it may be convenient to define thegrid shift s= (s1,s2,s3)⊺ in the basis (a∗m,b∗m,c∗m). Thereforewe have s̃= Ps.The indexing of the grid points in the reciprocal primitivecell (ã∗p , b̃∗p , c̃∗p) is a trivial task. Indeed, in the phonopy andphono3py codes, each grid point is bijectively mapped to aninteger p byp= m1 + n1m2 + n1n2m3, mi ∈ {0,1, . . . ,ni − 1}. (66)With this definition, 0⩽ p< n1n2n3.The q points generated this way are located within thereciprocal cell (ã∗p , b̃∗p , c̃∗p). They could be shifted to the recip-rocal cell (a∗p ,b∗p ,c∗p) by reciprocal lattice vector translationsif needed. Notice however that to obtain the integer p through9J. Phys.: Condens. Matter 35 (2023) 353001 Topical Reviewequation (66) for a general integer triplet (m1,m2,m3)⊺, mod-ulo n= (n1,n2,n3)⊺ is required to locate the point within thecell (ã∗p , b̃∗p , c̃∗p). Indeed, (m1,m2,m3)⊺ and (m1 +G1n1,m2 +G2n2,m3 +G3n3)⊺ indicate different locations in q space,although they are equivalent points due to periodicity.6.4. Symmetry of generalized regular gridsFor a Γ centred grid (si = 0), the coordinates of a q point inthe basis (a∗p ,b∗p ,c∗p) are given byq= QD−1m+G, (67)where m= (m1,m2,m3)⊺ and the reciprocal lattice vector Gis chosen to bring qi in the interval [0,1). If the regular gridis a traditional one, we simply have Q= P= 1. As shown byequation (18), the image of this q point through an operationof the crystallographic point group is given byq ′ = S−⊺QD−1m+S−⊺G+G ′ (68)where we have added a reciprocal lattice vector G ′ to shiftq ′i of the image point in the interval [0,1). If q ′ belongs tothe grid we have defined, for all S−⊺ in the crystallographicpoint group, we will say that the grid is invariant under thecrystallographic point group. If it is so, q ′ can be written asq ′ = QD−1m ′ +G ′ ′ (69)with m ′i ∈ {0,1, . . . ,ni − 1}. Comparing equations (68)and (69), we obtainm ′ = (QD−1)−1S−⊺(QD−1)m (70)+(QD−1)−1(S−⊺G+G ′ −G ′ ′) (71)Because Q and S−⊺ are unimodular, and D contains the num-ber of divisions along (ã∗p , b̃∗p , c̃∗p), the last term is always areciprocal lattice vector. The matrixS̃−⊺= (QD−1)−1S−⊺(QD−1) (72)is the matrix representation of S in the (ã∗m, b̃∗m, c̃∗m) basis. Itsdeterminant is always 1 or −1, but for m ′ to be integer, it hasto have integer entries, and therefore be unimodular. This ischecked in the phono3py code, and if it is true for all S−⊺ inthe crystallographic point group, the regular grid follows thecrystallographic point group, and we consider it is properlydefined. The last equation can also be written asP−1S̃−⊺P=MgS−⊺M−1g (73)and an equivalent strategy is to check that both sides areunimodular.Introducing the subgrid shift s (see section 6.3) can furtherbreak the symmetry of the regular grid. In equation (67) mis replaced by m+Ps and m ′ by m ′ +Ps in equation (69).Requiring the shift to be invariant (mod 1) under an operationof the crystallographic point group, S̃−⊺Ps= Ps (mod 1), weobtain the conditionP−1S̃−⊺Ps=MgS−⊺M−1g s= s (mod 1) (74)against all S−⊺ in the crystallographic point group.6.5. Double grid for subgrid shiftTo satisfy the crystallographic symmetry, the subgrid shift si isnormally chosen to be either 0 or 1/2. It is better to treat gridpoint arithmetic by integers for the computer implementationand its performance. This is realized by doublingm and s. Thiswell-known technique, e.g. presented in [3], is directly usablefor the generalized regular grid in section 6.2. As a choice, wedefine it bymd = 2(m+Ps), (mod 2n), (75)and the q points are given like in equation (67) byq= QD−1md/2+G. (76)The symmetry operation is implemented asmd ′ = S̃−⊺md (mod 2n). (77)The index of the grid point md ′ is obtained by equation (66)after recovering m ′ using equation (75),m ′ = (md ′ − 2Ps)/2 (mod )n. (78)In the phono3py code, an integer vector of 2s is used forthe implementation, instead of s, to avoid floating pointarithmetic.6.6. Grid points in BZIn this section, the BZ grid points are defined to include dif-ferent q points on the BZ surface that may be equivalent gridpoints by reciprocal lattice translations, in addition to the gridpoints inside the BZ.For a Γ centred grid, in the (ã∗m, b̃∗m, c̃∗m) basis, each BZ gridpoint is represented by the integer tripletmBZ =mBZ1mBZ2mBZ3=m1 +Gm1 n1m2 +Gm2 n2m3 +Gm3 n3 , Gmi ∈ Z, (79)where Gm is chosen to minimize the norm of q=(ã∗p , b̃∗p , c̃∗p)D−1mBZ = (ã∗m, b̃∗m, c̃∗m)mBZ. Therefore, mBZ iswritten asmBZ =m+D argminGm(|(ã∗p , b̃∗p , c̃∗p)(D−1m+Gm)|). (80)Multiple mBZ can be found when the grid point is on the BZsurface.In the phono3py code, equation (80) is implemented as fol-lows. As the first step, reduced basis vectors of the reciprocal10J. Phys.: Condens. Matter 35 (2023) 353001 Topical Reviewprimitive cell are obtained by using the Niggli reduction or anyreduction scheme. This is written using the change-of-basismatrixM∗p→r as(a∗r ,b∗r ,c∗r ) = (a∗p ,b∗p ,c∗p)M∗p→r, (81)whereM∗p→r is unimodular. This givesq= (a∗r ,b∗r ,c∗r )(M∗p→r)−1QD−1[m+DGm] (82)= (a∗r ,b∗r ,c∗r )[(M∗p→r)−1QD−1m+G]. (83)The second line defines the reciprocal lattice vector G. G isdivided into two pieces,G= Ḡ+∆G. (84)Nearest integers of−(M∗p→r)−1QD−1m are stored in Ḡ, whichis formally written asḠ=−nint[(M∗p→r)−1QD−1m] (85)to bring the following q̄ closer to the origin,q̄= (a∗r ,b∗r ,c∗r )[(M∗p→r)−1QD−1m+ Ḡ]. (86)∆G that minimizes |q̄+(a∗r ,b∗r ,c∗r )∆G| is searched in∆Gi ∈{−2,−1,0,1,2}. Finally, mBZ is obtained asmBZ =m+DQ−1M∗p→r(Ḡ+∆G). (87)As written above, m and p are in one-to-one correspond-ence and multiple mBZ can be found for each m or p. In thephono3py code, {{mBZ(p)}|0⩽ p< det(D)} is calculated andstored at the initial step.With the subgrid shift, m in the equations above arereplaced by m+Ps. Either with or without the subgrid shift,mBZ is determined using the double grid technique presentedin section 6.5. From stored mBZ, the BZ q points are obtainedbyq= (a∗p ,b∗p ,c∗p)QD−1(mBZ +Ps). (88)6.7. Irreducible grid pointsUsing symmetry properties of phonons, e.g. equation (25),the computation required for the phonon calculations can bereduced. Indeed, many phonon properties result from the BZintegration of functions of a single q variable. Collecting thevalues of those functions for the q points in the irreduciblepart of the BZ only may speed up the calculation and also savememory space.When the BZ integration is performed on a reg-ular grid, the irreducible q points are obtained from the irre-ducible grid points. In this section, we explain how to obtainthose irreducible grid points, where it is assumed that the reg-ular grid satisfies the symmetry as described in section 6.4.A 2D example of BZ is presented in figure 8. A symmetryoperation acts on a q point as q ′ = Sq. The star of q is the setFigure 8. (a) 2D Brillouin zone (BZ). Shaded area is a choice of theirreducible BZ. (b) 6× 6 regular grid on the BZ including the BZsurface. Filled and open circle symbols show the grid pointsbelonging to the irreducible BZ under the C4v symmetry and theother grid points, respectively. (c) The irreducible grid points maynot be located in a connected space.of those points written as {Sq|S ∈ P}, where one q point rep-resents its star. The irreducible BZ is the set of representativeq points of all stars in the BZ. In figure 8(a), the 2D BZ has theeight-fold symmetry of C4v, and the irreducible BZ is depictedby the shaded area.In a phonon calculation, q points are sampled on a regu-lar grid, and the irreducible grid points are defined as the gridpoints that belong to the irreducible BZ as shown in figure 8(b).In practical calculations, the irreducible grid points consideredmay however not be located in a connected space, as shown infigure 8(c). This simply means that the selected representativepoint in the star is located somewhere else.In the phonopy and phono3py codes, the irreducible gridpoints are obtained as follows. Each grid point indexed by p isrepresented in the double grid of equation (75). All symmetryoperations S̃−⊺of the crystallographic point group are appliedto the grid point as given by equation (77). The grid pointindex of the rotated grid point, p(S), is recovered success-ively applying equations (78) and (66). The minimum valuein {p(S)|p(S)⩽ p,S ∈ P} is chosen as the irreducible gridpoint. The mappings of the grid point indices p{S}−−→ p arestored for all grid points in 0⩽ p< det(D), and the uniqueelements of {p} give the irreducible grid points of the regu-lar grid. It is important to perform these operations only byintegers for computational efficiency.6.8. q point tripletsIn the phono3py code, triplets of q points, (q,q ′,q ′ ′), are con-sidered for the computation of the phonon–phonon interaction.This interaction is therefore not a function of a single q vari-able. The number of grid points det(D) becomes large whendense sampling is used, e.g. det(D) = 3003 = 27× 106 in thestudy of [20]. Then, the number of triplets of q points becomes3009 ∼ 1022, which is a really huge number for numericalcomputations and also storing in the memory of computers.Due to lattice translational symmetry, elements of the threephonon interaction strength can be non-zero only if q+q ′ +q ′ ′ =G [32, 33]. This is used to reduce the computationof the interactions for all triplets (q,q ′,q ′ ′) to all pairs of11J. Phys.: Condens. Matter 35 (2023) 353001 Topical Reviewpoints (q,q ′). To satisfy this condition, the Γ centred regu-lar grid is used in the phonon–phonon interaction calculationof the phono3py code. Moreover, the crystallographic pointgroup operations, time-reversal symmetry, and permutationsymmetry for each pair of q and q ′ are used to reduce fur-ther the number of interactions to be computed. However, evendoing so, the number of the combination of two q points is stilllarge. Therefore, in this section, we describe the strategy usedto make those computations practical.As mentioned in section 6.7, many phonon related proper-ties are written as sum of functions of a single variable q, suchas∑q f(q). For example, under the relaxation time approxim-ation, the LTC computed in the phono3py code can be writtenasκ=1NVc∑qντqνCqνvqν ⊗ vqν , (89)where Cqν and vqν denote the phonon mode heat capacity andgroup velocity, respectively,Cqν = kB(ℏωqνkBT)2 exp(ℏωqν/kBT)[exp(ℏωqν/kBT)− 1]2, (90)vqν =∂ωqν∂q. (91)In equation (90), ℏ and kB denote the reduced Planck constantand the Boltzmann constant, respectively, and T is the tem-perature. τqν in equation (89) is the phonon lifetime, and thereciprocal of τqν is calculated from the phonon–phonon inter-action strength in the phono3py code. The explicit equationused for τqν = 1/2Γqν(ωqν) is given at equation (111).The LTC of equation (89) is conveniently computed iterat-ing over irreducible q points. Indeed, the mode heat capacityand the lifetime have the same symmetry as the phonon bandstructure, equation (25),CSqν = Cqν , (92)τSqν = τqν , (93)and the derivative of equation (25) givesvSqν = Svqν . (94)Therefore, the summation in equation (89) can be reduced toa summation over the irreducible grid points. Denoting q theirreducible points, we obtainκ=1NVc∑qντqνCqν(∑SSvqν ⊗Svqν)m(q)|{S}|(95)where the last factor is the number of branches in the star of q,m(q), divided by the cardinality of the crystallographic pointgroup, |{S}|.When each calculation of τqν is computationally demand-ing, the computation of τqν at different q points may be dis-tributed over multiple or many computer nodes. Therefore, itFigure 9. A strategy to choose a q point triplet. The same 2DBrillouin zone (BZ) as figure 8 is used. The star symbols depictthree q points in the BZ under the constraint of q+ q ′ + q ′ ′ = G.(a) q is chosen from the irreducible grid points of C4v symmetry. qworks as the stabilizer of the subgroup of C4v symmetry. (b) q ′ ischosen in the irreducible grid points of the subgroup. (c) q ′ ′ = G−q− q ′ where G of the shortest |(a∗,b∗)G|, i.e. G= (0,0) in thisexample, is chosen.is convenient to have a set of q point triplets at fixed q in{(q,q ′,q ′ ′)}. As shown in equation (95), in this strategy, q ischosen from the irreducible BZ. Now, because q must be keptfixed, it is not the full set of operations of the crystallographicpoint group which should be used to reduce the summationover q ′, but a subgroup which lets q invariant. This subgroupis known as the point group of q,Pq = {S|Sq= q (mod G), S ∈ P}. (96)Consequently, q ′ is chosen in the irreducible part of the BZdefined by Pq. Finally, q ′ ′ is computed as q ′ ′ =G−q−q ′whereG is a reciprocal lattice vector chosen to shift q ′ ′ withinthe BZ of the same origin. For phonons, one of (q,q ′,q ′ ′) and(q,q ′ ′,q ′) is chosen due to the permutation symmetry.In the phono3py code, q, q ′, and q ′ ′ are taken from theBZ of the same origin. When some of them are on the BZsurface, they are chosen among their translationally equivalentpoints on the BZ surface so that the triplet minimizes |G| ofq+q ′ +q ′ ′ =G.This triplet search is implemented using the irreducible gridpoints and the BZ grid points described in sections 6.7 and 6.6,respectively. The q points q, q ′, and q ′ ′ are represented by theBZ grid pointsmBZq ,mBZq ′ , andmBZq ′ ′ , respectively. As shown infigure 9(a), mBZq is chosen in the irreducible grid points. Thepoint mBZq breaks the crystallographic point-group P. UsingS̃−⊺given by equation (72), the point group of q is obtainedasPq ={S|mBZq = S̃−⊺mBZq (mod n), S ∈ P}. (97)As shown in figure 9(b), mBZq ′ is sampled from the irreduciblegrid points of Pq. The third grid point is given by mBZq ′ ′ =DG−mBZq −mBZq ′ . Finally, mBZq , mBZq ′ , or mBZq ′ ′ may be shif-ted to minimize |G| if some of them are on the BZ surfaceas explained for (q,q ′,q ′ ′) above.12J. Phys.: Condens. Matter 35 (2023) 353001 Topical ReviewFigure 10. Implementation of the linear tetrahedron method in thephonopy and phono3py codes. This routine is provided as a modulethat returns integration weights. In the module, each layer dependson the lower layers. Crystal symmetry is treated outside of thismodule.Figure 11. A scheme to divide a parallelepiped microzone into sixtetrahedra with the same volume [3]. The eight vertices of theparallelepiped are shared by 6, 6, 2, 2, 2, 2, 2, and 2 tetrahedra,respectively.7. Tetrahedron methodOnce implemented, tetrahedron methods are easy to use androbust. Moreover, when the code is written in a modularway, it can be reusable for different kinds of BZ integrations.The phonopy and phono3py codes employ a linear tetrahed-ron method where several techniques are picked and mixedfrom the various reports, notably, by MacDonald et al [2] andBlöchl et al [3]The implementation in the phonopy and phono3py codes isstructured as schematically illustrated in figure 10. The bottomof the structure is the regular grid. The parallelepiped formedby the basis vectors of the reciprocal primitive cell is dividedinto many smaller parallelepiped microzones on the regulargrid. Everymicrozone has the same shape, defined by the basisvectors (a∗m,b∗m,c∗m), as presented in section 6. In the case ofthe traditional regular grid, equation (56) defines the micro-zones to be considered, while for the generalized regular gridequation (59) is used. Each of the microzones is divided intosix tetrahedra in the same way as shown in figure 11.The integral over the BZ to be performed is then regardedas a sum of contributions from those tetrahedra. The functionto be integrated is linearly interpolated within each tetrahed-ron (see section 7.2) and therefore the integration within eachtetrahedron can be done analytically. To interpolate linearly atridimensional function one needs four values. Each tetrahed-ron has four vertices. Therefore the values of the function tobe integrated on the vertices of the tetrahedra are used to buildthe interpolations.Different recipes can be used to implement this program.We follow [2] as well as [3] to obtain integration weights forgrid points rather than for the tetrahedra themselves. Indeed,as it will be shown, it is possible to rearrange the contributionsfrom tetrahedra to contributions from the grid points. This isuseful to make the tetrahedron method easy use, and to makeit similar to smearing methods. In addition, this way, the sym-metry of the regular grid described in section 6.4 is appliedstraightforwardly to the integration weights.7.1. Division into six tetrahedraThe scheme to divide themicrozone into six tetrahedra follows[3] reported by Blöchl et al We choose a shortest main diag-onal of the parallelepiped (a∗m,b∗m,c∗m). Then the six tetrahedraare selected sharing this main diagonal as shown in figure 11.7.2. FunctionsIn this section we introduce the functions g(ω), n(ω), I(ω),and J(ω) to be computed using the linear tetrahedron method.The notation roughly follows [2]. Moreover, the band index νis not written explicitly, since each band can be treated inde-pendently. It is therefore understood that a sum over the bandsshould be performed at the end of the calculation.The density of states g(ω) is written asg(ω)≡ˆBZd3qδ(ω−ωq) (98)≈ Vt6N∑i=1g(ω,ωi1,ωi2,ωi3,ωi4)≡ Vt6N∑i=1gi. (99)The definition is written at the first line, the approximationobtained from the linear tetrahedron method at the second lineand an obvious definition of gi, the contribution of tetrahedroni, at the third line. Here i is an index running throughout the6N tetrahedra, and Vt is the volume of a tetrahedron. The val-ues of frequency at the vertices are assumed to be arranged inascending order, ωi1 < ωi2 < ωi3 < ωi4.The integrated density of states, or number of states func-tion, n(ω), is written following the same conventions,n(ω)≡ˆ ω−∞dω ′g(ω ′) (100)≈ Vt6N∑i=1n(ω,ωi1,ωi2,ωi3,ωi4)≡ Vt6N∑i=1ni. (101)13J. Phys.: Condens. Matter 35 (2023) 353001 Topical ReviewWeighted density of states frequently appears in phononcalculations. They are defined byI(ω) =ˆBZd3qF(q)δ(ω−ωq)≈ Vt6N∑i=1gi4∑k=1Ik(ω,ωi1,ωi2,ωi3,ωi4)Fik≡ Vt6N∑i=1gi4∑k=1I ikFik, (102)where F(q) is a function of q and we used the notationFik = F(q ik) at the vertex k of the tetrahedron labelled by i.I ik = Ik(ω,ωi1,ωi2,ωi3,ωi4) are given in appendix A. We noticein this equation an additional summation over the vertices ofeach tetrahedron. This comes from the functionF, which is lin-early interpolated within each tetrahedron. When F= 1, thisequation reduces to the density of states.Finally, the integral of I(ω) over frequency, J(ω), is writtenasJ(ω) =ˆ ω−∞dω ′I(ω ′) =ˆBZd3qF(q)θ(ω−ωq)≈ Vt6N∑i=1ni4∑k=1Jk(ω,ωi1,ωi2,ωi3,ωi4)Fik≡ Vt6N∑i=1ni4∑k=1J ikFik, (103)where J ik = Jk(ω,ωi1,ωi2,ωi3,ωi4) are given in appendix A.In appendix A, the formulae for the coefficients appearingfor the quantities, gi, ni, I ik, Jik, are stated from [2]. To derivethose formulae, the calculation of J(ω) is considered first, andI(ω) is obtained through differentiation. n(ω) and g(ω) arejust special cases for F= 1. The θ(ω−ωq) function in J(ω)shows that computing the contribution from one tetrahedronis related to the estimation of the volume ωq < ω, which canbe described from an intersection of a tetrahedron with a planeωq = ω, since ωq is now assumed to be a linear function of q.As shown in figure 12 the volume will assume different shapesdepending on the value of ω. Different formulae are thereforeobtained for the different cases to be considered.7.3. Integration weightsIn equation (102), I(ω) is expressed as a summation over 6Ntetrahedra and four vertices. Following Blöchl et al in [3], thiscan be rearranged as a sum over grid points. This scheme hidesthe cumbersome handling of tetrahedra from outside the lin-ear tetrahedron method module as shown in figure 10. Onceintegration weights, that are described in this section, are com-puted, the data structure is shared with smearing methods andthe symmetry handling is made easy.Figure 12. Different linear cuts by ω plane of a tetrahedron.(a) ωi1 < ω < ωi2 < ωi3 < ωi4, (b) ωi1 < ωi2 < ω < ωi3 < ωi4, and(c) ωi1 < ωi2 < ωi3 < ω < ωi4 [2].I(ω) can be rewritten asI(ω)≈ VtN−1∑p=024∑ip=1gip4∑k=1Iipk Fipk δ(ip,k),p, (104)where p denotes the grid points (see section 6.3), and ip is theindex of a tetrahedron in the 24 tetrahedra that are sharing thegrid point p (see figure 11 for 24= 6+ 6+ 2+ 2+ 2+ 2+2+ 2). The composite index (ip,k) indicates a grid point, andδ(ip,k),p selects the vertex which is identical to the grid point pin the four vertices of the tetrahedron ip.Since Fipk δ(ip,k),p is only non-zero when (ip,k) is the pointp, we write Fipk δ(ip,k),p = Fpδ(ip,k),p. The summation (104)becomesI(ω)≈ VtN−1∑p=0Fpwp, (105)wherewp =24∑ip=1gip4∑k=1Iipk δ(ip,k),p. (106)I(ω) has therefore been expressed as a summation over gridpoints. This integral is given by the weighted sum of the valuesof the function F at points p, Fp, with weights wp given byequation (106).A 2D illustration of the summation (105) is presented infigure 13. Figure 13(a) shows that each parallelogram is cutin two triangles, and the four vertices of each parallelogramare shared by 2, 2, 1, and 1 of those triangles. Six trianglesshare the grid point marked by the circle symbol as shown infigure 13(b). This shows that only six terms are summed up inthe summation (106) because of δ(ip,k),p. For the tetrahedra in3D, 24 terms are summed up at each grid point in the summa-tion (106).To compute equation (106), it is convenient to representthe positions of grid points using shifts from the grid point p,mp+∆m. The 2D example is depicted in figure 14. The maindiagonal (0,0)–(1,1) is chosen. For this main diagonal, the sixtriangles that share the central point are uniquely determined.14J. Phys.: Condens. Matter 35 (2023) 353001 Topical ReviewFigure 13. 2D schematic illustration of rearrangement ofsummation (104). (a) Each parallelogram is cut into two triangles.Four vertices of the parallelogram are shared by 2, 2, 1, and 1triangles, respectively. (b) Six triangles (shaded) that share a gridpoint (circle symbol) contribute to the integration weight of the gridpoint.Figure 14. 2D schematic illustration of grid point shifts∆m fromthe focused grid point (centre) to the neighbours.They are represented by the following ∆m of the apices:(0,0)–(0,1)–(1,1),(0,0)–(1,1)–(1,0), . . ..For the linear tetrahedronmethod, there are four choices forthe main diagonal, and for each choice of the main diagonal,the ∆m of the vertices of the 24 tetrahedra shared by a gridpoint are predetermined. Therefore, the datasets of the 24 tetra-hedra of the four main diagonals are predetermined and hard-corded in the phonopy and phono3py codes as 4× 24× 4× 3integer array whose elements are in {−1,0,1} for the regulargrid.For the generalized regular grid, the shifts −1 or 1 cor-respond to the coordinates of the neighbouring points in the(a∗m,b∗m,c∗m) basis. The coordinates of those shifts are howevereasily obtained in the (ã∗m, b̃∗m, c̃∗m) basis used in the calcula-tion, ∆mG = P∆m, since(a∗m,b∗m,c∗m)∆m= (a∗m,b∗m,c∗m)P−1P∆m= (ã∗m, b̃∗m, c̃∗m)∆mG, (107)where equation (63) is used in the last equation.7.4. Use of symmetry of grid pointsThe symmetry of grid points can be applied to the lineartetrahedron method. The symmetry is used to prepare theinput dataset of a grid point required by the linear tetrahed-ron method module as illustrated in figure 10. Once the inputdataset is made, symmetry information is unnecessary for thelinear tetrahedron method module.In this section, we denote ωp ≡ ωik and the mapping ofthe grid point p{S}−−→ p defined in section 6.7 is written asp= f(p). Under the condition ωp = ωf(p), which is satisfiedby the phonon frequencies, the integration weights fulfil therelationwp ≈ wf(p). (108)This is an approximation since the set of 24 tetrahedra maynot follow the crystallographic point group, i.e. the rotatedtetrahedron may not be mapped to any of the 24 tetrahedra.However, equation (108) is a very good approximation.When Fp has the symmetry Fp = Ff(p), equation (105) canbe written asI(ω)≈ VtN−1∑p=0Fpwp ≈ VtN−1∑p=0Ff(p)wf(p). (109)The right hand side of Approx. (109) can therefore be com-puted from the values at irreducible grid points only. Definingthe multiplicity of the irreducible points as m(p) = |{p|p=f(p), ∀p}|, the right hand side of Approx. (109) is written asI(ω)≈ Vt∑{p}Fpwpm(p). (110)Herem(p) is equivalent tom(q) in equation (95), i.e. the num-ber of branches in the star of q.7.5. Triplets integration weightsAs shown in [33], the phonon linewidth computed by thephono3py code, due to phonon–phonon interactions, is givenbyΓqν(ω) =18πℏ2∑q ′q ′ ′∑ν ′ν ′ ′∣∣Φ−qν,q ′ν ′,q ′ ′ν ′ ′∣∣2{(nq ′ν ′ + nq ′ ′ν ′ ′ + 1)δ(ω−ωq ′ν ′ −ωq ′ ′ν ′ ′) (111)+(nq ′ν ′ − nq ′ ′ν ′ ′) [δ(ω+ωq ′ν ′ −ωq ′ ′ν ′ ′)− δ(ω−ωq ′ν ′ +ωq ′ ′ν ′ ′)]} . (112)15J. Phys.: Condens. Matter 35 (2023) 353001 Topical ReviewIf we consider q to be fixed, the summations to be per-formed over the BZ have the formKq(ω) =1N∑q ′q ′ ′Fq,q ′,q ′ ′∆(q+q ′ +q ′ ′)δ(ω−ωq ′,q ′ ′)(113)where Fq,q ′,q ′ ′ is an arbitrary function of q, q ′ and q ′ ′.ωq ′,q ′ ′is a function of q ′ and q ′ ′. For example, in equation (111), tocompute the contribution from the first delta function we needto consider, ωq ′,q ′ ′ = ωq ′ +ωq ′ ′ .In the above equation, ∆(q+q ′ +q ′ ′) means 1 whenq+q ′ +q ′ ′ =G otherwise 0. In equation (111) this factor isincluded inΦqν,q ′ν ′,q ′ ′ν ′ ′ [33]. For given q, q ′, and q ′ ′ in theBZ, q ′ ′ is chosen so that |G| is smallest. At fixed q, for givenq ′, only a single q ′ ′ contributes. This allows to eliminate thesummation over q ′ ′, therefore we obtainKq(ω) =1N∑q ′Fq,q ′,G−q−q ′δ(ω−ωq ′,G−q−q ′) (114)=Vc(2π)3ˆBZd3q ′Fq,q ′,G−q−q ′δ(ω−ωq ′,G−q−q ′).(115)This quantity has therefore the form of the function I(ω)we considered previously. To apply the linear tetrahedronmethod to Kq(ω) we need to know the value of the functionsFq,q ′,G−q−q ′ and ωq ′,G−q−q ′ at the vertices q ′ +∆q ′ aroundthe point q ′. They are given byFq,q ′+∆q ′,G−q−q ′−∆q ′ = Fq,q ′+∆q ′,q ′ ′−∆q ′ , (116)ωq ′+∆q ′,G−q−q ′−∆q ′ = ωq ′+∆q ′,q ′ ′−∆q ′ . (117)It means that∆q ′ ′ =−∆q ′ for the corresponding neighbour-ing point of q ′ ′. Using the valuesωq ′+∆q ′,q ′ ′−∆q ′ as the inputdataset of the linear tetrahedron method, the triplets integ-ration weights are computed in the same way as written insection 7.3.8. Random displacement generationOnce phonon frequency ωqν and eigenvector Wκα(qν) areobtained in a supercell, they can be used to generate atomicconfigurations relevant at finite temperatures. This can beachieved using equation (27) and writingQ(qν) =1√2(QR(qν)+ iQI(qν)), (118)where both QR(qν) and QI(qν) are real variables. It is shownin [34] that they are real normal coordinates, fulfilling theequation of the harmonic oscillator. Their probability densityis therefore a Gaussian,P[QR/I(qν)] =1√2πσ2qνexp[−QR/I(qν)22σ2qν], (119)withσ2qν =ℏ2ωqνcothℏωqν2kBT=ℏ2ωqν(1+ 2nqν), (120)where nqν = (eℏωqνkBT − 1)−1 is the Bose–Einstein distributionfunction.Substituting equation (118) into equation (27) and remem-bering Q(−qν) = Q∗(qν), we obtainulκα =1√Nmκ{ ∑q∈A,νQ(qν)Re[ulκα(qν)]+1√2∑q∈B,ν[(QR(qν)+ iQI(qν))ulκα(qν)+ (QR(qν)− iQI(qν))u∗lκα(qν)]}(121)=1√Nmκ{ ∑q∈A,νQ(qν)Re[ulκα(qν)]+√2∑q∈B,ν[QR(qν)Re[ulκα(qν)]−QI(qν)Im[ulκα(qν)]]}. (122)Here, A is a set of q points commensurate with the supercell,where Q(qν) becomes real (e.g. Γ point) or q=−q+G, andB includes other commensurate q points in the one side of theBZ and q ̸=−q+G. Namely, only one of q or −q should beincluded in the summation over q ∈ B [34].Finally, generating random numbers followingequation (119) and using them for QR(qν) and QI(qν) inequation (122), random displacements {ulκα} in a supercellcan be generated systematically at a given temperature.The random structures generated in this way can be used toinvestigate the impact of phonon excitation on various phys-ical properties, such as electronic structure and magnetism, atfinite temperatures. Indeed, the investigated quantity is com-puted for each random structure generated, and the averageof values obtained for those random structures is computed tomimic the effect of temperature.We note that this approach is physically valid only whenall phonon modes are dynamically stable (ω2qν > 0) at allcommensurate q points. This condition is usually satisfiedin many compounds where lattice anharmonicity is not verystrong. However, in some anharmonic systems, includinghigh-temperature metastable phases, the harmonic approxim-ation may yield unstable phonon modes.Even when all phonons are stable, their frequencies maydeviate from experimental values significantly. In such cases,it is recommended to use anharmonic phonon frequenciesand eigenvectors at finite temperatures for equations (120)and (122). The anharmonic frequencies and eigenvectorsrenormalized at finite temperatures can be obtained, forexample, by performing a self-consistent phonon calculation16J. Phys.: Condens. Matter 35 (2023) 353001 Topical Reviewusing either the real-space stochastic approaches, [35–40]which are often termed the stochastic self-consistent harmonicapproximation (SSCHA), or themomentum-space implement-ation which uses anharmonic force constants [41, 42]. Anapplication of SSCHA performed with random displacementsgenerated by the phonopy code coupled with the ALM forceconstants calculation code [43] is found in [20].9. Phonon band unfoldingTo draw a phonon band structure for a system with defectscomputed in a supercell approach, band unfolding techniquemay be useful. Since even one point defect breaks the peri-odicity of a crystal, a certain approximation is necessary toachieve this representation. In this section, an implementationof the method proposed by Allen et al [4] is explained.In figure 15, a 2D schematic illustration of a supercell witha vacancy is presented. It is analyzed by presuming the cor-responding perfect supercell with the same lattice vectors.Equation (30) gives the relationship between the basis vectorsof the primitive cell and those of the supercell, which indicatesthat the perfect supercell contains N= det(Mp→s) primitivecells. The reciprocal basis vectors are given by equation (46).The BZ volume of the primitive cell is N times larger than thatof the supercell as shown in figure 16.For this analysis to be possible, we assume there exists aone-to-one correspondence between the atoms and vacancysites in the defective supercell and the atoms in the perfectsupercell. The atoms and vacancies in the defective supercellcan therefore be labelled by the composite index lκ in the sameway as in the perfect supercell. It means that vacancy sites inthe defective supercell must have respective atoms in the per-fect supercell. The atoms in the defective supercell can nev-ertheless have slightly different positions than in the perfectsupercell. The 2D schematic illustration we obtain is shown infigure 15.We write the coordinates of lattice points in the supercellasRl = (ap,bp,cp)l= (as,bs,cs)̃l, (123)where l are integers, while l̃ may not be.Let T̂j be the jth lattice translation operator which is definedon a function as followsT̂jf(R) = f(R−Rlj). (124)In this section, l correspond to lattice vectors within the per-fect supercell. There are N such lattice vectors and thereforeN lattice translation operations. T̂j maps a lattice vector withindex li in the perfect supercell to a lattice vector with indexlk. The mapping between indices is determined froml̃iT̂j−→ l̃k = l̃i + l̃j (mod 1), (125)where (mod 1)means that the lattice points translated outsidethe supercell are brought back inside the supercell by a super-cell lattice translation. Themapping obtained this way is storedFigure 15. 2D schematic illustration of a supercell with a vacancy(left) and the corresponding perfect supercell (right). The atoms aredepicted by open circle symbols. The vacancy is marked by “V”.The perfect supercell contains 2× 2 primitive cells whose latticepoints (l1, l2, l3, l4) are indicated by filled circle symbols. Positionsof all atoms in the defective supercell can deviate slightly from theequilibrium positions of the perfect supercell. At unfolding, thephase shift of q⋆ ·Rlj is multiplied to the phonon eigenvector of thedefective supercell at each primitive lattice translation as given inequation (135).Figure 16. Unfolding of 2D Brillouin zone (BZ) of the supercell(q⋆1 ) to 2× 2 BZs of the supercell (q⋆1 ,q⋆2 ,q⋆3 ,q⋆4 ). Gray levels ofshaded backgrounds indicate BZs of the supercell which belong todifferent q⋆ points.in an N×N permutation matrix. In figure 15, the 2D supercellmodel is made of 2× 2 primitive cells. For example, the trans-lation by ap brings l1 to l2. One more translation by ap bringsl2 to l1 due to the periodicity of the supercell. In the following,the mapping liT̂j−→ lk is written compactly as lk = Tj(li).The first step of this method is to compute the phononeigenvectors of the defective supercell. The Bloch wave of aneigenvector given by equations (19) and (20) is written asΨlκ(q̃ν̃) = eiq̃·RlκWlκ(q̃ν̃), (126)where Wlκ(q̃ν̃) is the phonon eigenvector of the defectivesupercell, with band index ν̃, and at the wave vector q̃ in theBZ of the supercell. Eigenvector elements corresponding tovacancies are treated as zeros.We consider that T̂j is applied to the Bloch wave asT̂jΨlκ(q̃ν̃)≈ΨT−1j (l)κ(q̃ν̃)= eiq̃·(Rlκ−Rlj )WT−1j (l)κ(q̃ν̃). (127)17J. Phys.: Condens. Matter 35 (2023) 353001 Topical ReviewThe first line of equation (127) may not be consistent withequation (124) since positions of atoms in the defectivesupercell can be displaced slightly from those in the perfectsupercell.The next step is to unfold the q̃ point from the BZ of thesupercell to theN q points in the BZ of the primitive cell whichare equivalent up to vectors of the reciprocal supercell lattice.This is written asq= q⋆ + q̃, (128)where q⋆ is an integer linear combination of (a∗s ,b∗s ,c∗s ). Theschematic illustration for the BZ of the 2D supercell is shownin figure 16. The q̃ points in the BZ of the supercell are shiftedby q⋆ inside the BZ of the primitive cell.Applying periodic boundary conditions on the supercell,the translation operations {T̂j} form an Abelian group of orderN with irreducible representations e−iq⋆·Rlj . The great ortho-gonality theorem then gives the relation1NN∑j=1eiq⋆·Rlj = δq⋆,0. (129)This result is also proved algebraically in appendix B, togetherwith the dual relation1NN∑i=1e−iq⋆i ·(Rlj−Rlj ′)= δj,j ′ . (130)where q⋆i are the N vectors q⋆ within the BZ of the primitivecell. We define the operatorP̂(q̃,q⋆) =1NN∑j=1ei(q̃+q⋆)·Rlj T̂j. (131)It can easily be shown that it is a projection operator, whichselects the part of Ψlκ(q̃ν̃) that is also a Bloch function ofthe primitive lattice with wave vector q= q⋆ + q̃. Indeed, inappendix C we obtainT̂i P̂(q̃,q⋆)Ψlκ(q̃ν̃) = e−i(q̃+q⋆)·Rli P̂(q̃,q⋆)Ψlκ(q̃ν̃). (132)Moreover, because from equation (130),∑Ni=1 P̂(q̃,q⋆i ) = 1,the decomposition of Ψlκ(q̃ν̃) into the N functionsP̂(q̃,q⋆i )Ψlκ(q̃ν̃), i = 1, . . . ,N is unique,Ψlκ(q̃ν̃) =N∑i=1[P̂(q̃,q⋆i )Ψlκ(q̃ν̃)]. (133)The projection by equation (131) and the decompositionby equation (133) indicate that the weight of state havingthe Bloch symmetry q= q⋆ + q̃ of the primitive lattice in thesupercell Bloch state Ψliκ(q̃ν̃) is therefore defined asFigure 17. Phonon band unfolding of Al with a vacancy. Color scaleindicates unfolding weights of phonon modes depicted by filledcircle symbols, where those of the degenerated modes are summed.Phonon modes calculated for the defective supercell are presentedby + (plus) symbols. Unshaded and shaded backgrounds showBrillouin zones of the supercell that belong to q⋆ of Γ and L,respectively. Only a part of the weights are unfolded into these twocommensurate points among 32 commensurate points. Althoughthis unfolding technique requires no reference states, for thecomparison, the phonon band structure of the perfect supercell isshown by solid curves.w(q̃ν̃,q⋆) =∑iκ|P̂(q̃,q⋆)Ψliκ(q̃ν̃)|2=1N2∑iκ∣∣∣∣ N∑j=1WT−1j (li)κ(q̃ν̃)eiq⋆·Rlj∣∣∣∣2 (134)=1N∑iκN∑j=1W∗liκ(q̃ν)WT−1j (li)κ(q̃ν̃)eiq⋆·Rlj .(135)It can easily be checked that∑Ni=1w(q̃ν̃,q⋆i ) = 1 fromequation (130) and the unit normalization of the phonon eigen-vectors. Interestingly, this result does not require phononeigenvectors of the perfect crystal. The second equation isimplemented in the phonopy code though the third equationmay be more intuitive. WT−1j (li)κ(q̃ν̃) is the phonon eigen-vector elements belonging to the lattice point T−1j (li). Theunfolding weight indicates its correlation withWliκ(q̃ν) afterthe phase shift eiq⋆·Rlj . In the example of figure 15, the phaseshifts are integer multiples of π, since the 2D supercell modelis made of 2× 2 primitive cells.An example of the phonon band unfolding techniqueapplied to Al with a vacancy is shown in figure 17. Thecorresponding perfect supercell (32 atoms) is constructed by2× 2× 2 of the conventional unit cell (4 atoms). The primitivecell contains one atom and18J. Phys.: Condens. Matter 35 (2023) 353001 Topical ReviewMp→s =−2 2 22 −2 22 2 −2 . (136)Under fixed basis vectors of the supercell with a vacancy,internal atomic positions are relaxed. For the unfolding, qpoints are sampled along the BZ path L–Γ–L of the primitivecell (see figure 17). Those points can be uniquely decomposedas q= q̃+q⋆, where q̃ belongs to the BZ of the supercell, andq⋆ is a reciprocal lattice vector of the supercell. For the aboveexample, q⋆ is found to be Γ or L, depending on the positionof q along the path L–Γ–L. Those two possibilities are repres-ented using unshaded and shaded backgrounds in figure 17.The phonon modes of the defective supercell are calculatedat those q̃ points. The projection operator of equation (131)is q⋆ dependent and only a fraction of the phonon modes areunfolded along the L–Γ–L path, because the selected q⋆ pointsare only a subset of the commensurate q⋆ points in the summa-tion of equation (133). In figure 17, though we can see phononfrequencies are perturbed by the vacancy, the unfolded phononband structure roughly follows the phonon band structure ofthe perfect supercell.Substitutional defects can be managed straightforwardly.Interstitial defects can be handled in the same way asthe vacancy case by defining the corresponding supercelldifferently [44]. Suppose we have one interstitial defect in thesupercell model. It is treated as a usual atom in the primit-ive cell of the corresponding perfect supercell. In the defectivesupercell, the interstitial atom is located in one primitive cell,and in the other primitive cells, the interstitial sites are treatedas vacancies.10. ConclusionA scientific simulation code implements mathematical mod-els that are described with formulae. In those formulae, theinformation is often not detailed enough for direct computerimplementation, and the detailed knowledge can be implicit.Scientific software developersmay have to learn it from skilledpeople around, otherwise find it independently. To ease the lat-ter situation in the modern era of distributed software develop-ment, we aimed to provide details of the computational meth-ods employed in the phonopy and phono3py codes, keepingthe formalism as close as possible to the implementations.Not all of the computational methods described in thisreview are specific to the phonon calculation or the super-cell approach. For example, the handling of the regular gridsand the linear tetrahedron method presented are applicable toelectronic structure calculations. We therefore hope that thisreview will be a useful source of information for softwaredevelopers in condensed matter science.Data availability statementNo new data were created or analysed in this study.AcknowledgmentsThis work was supported by MEXT Japan through ESISM(Elements Strategy Initiative for Structural Materials) ofKyoto University and JSPS KAKENHI Grant NumbersJP21K04632 and JP21K03424. L C thanks the InstitutUniversitaire de France and the EXPLOR center hosted by theUniversity de Lorraine.Appendix A. Functions used in the lineartetrahedron methodIn this sectionwe give explicit formulae for gi, ni, I ik, and Jik [2].They assume different forms depending on the position of ωin the sequence ωi1 < ωi2 < ωi3 < ωi4. Also, they are expressedusing f inm and ∆inm which are defined asf inm ≡ (ω−ωim)/(ωin −ωim), finm+ f imn = 1, (A1)∆inm ≡ ωin −ωim. (A2)gi, ni, Ii, and Ji satisfy the following relations:dnidω= gi, (A3)d(ni J ik)dω= gi I ik. (A4)It may be useful to rememberf in ′m/∆inm = f inm/∆in ′m, (A5)df inmdω=1∆inm. (A6)A.1. ω < ωi1gi = 0. (A7)ni = 0. (A8)I ik = 0. (A9)J ik = 0. (A10)19J. Phys.: Condens. Matter 35 (2023) 353001 Topical ReviewA.2. ωi1 < ω < ωi2As shown in figure 12(a), the volume to consider is a tetra-hedron with vertices located at q1, qα = f i21q2 + f i12q1, qβ =f i31q3 + f i13q1, qγ = f i41q4 + f i14q1, andni = f i21 fi31 fi41. (A11)Inside this volume, F(q) is linearly interpolated, i.e. F(qα)≈f i21F(q2)+ f i12F(q1), . . ., and the contribution of ith tetrahed-ron to J(ω) is approximated asVtni [F(q1)+F(qα)+F(qβ)+F(qγ)]/4. By substituting it in equation (103), we obtainJi1 = (1+ f i12 + f i13 + f i14)/4. (A12)J ik = f ik1/4, k= 2,3,4. (A13)Using equations (A3) and (A4),gi = 3ni /(ω−ωi1) = 3f i21 fi31 fi41/(ω−ωi1)= 3f i21 fi31/∆i41. (A14)I i1 = ( f i12 + f i13 + f i14)/3. (A15)I ik = f ik1/3, k= 2,3,4. (A16)A.3. ωi2 < ω < ωi3In this case, the occupied part of the tetrahedron is the sumof the following three tetrahedra. Vertices of the first tetrahed-ron are q1, q2, qα = f i42q4 + f i24q2, qβ = f i32q3 + f i23q2, and itsvolume is Vi12αβ = Vt f i42 fi32. The vertices of the second tet-rahedron are q1, qα, qβ , qδ = f i41q4 + f i14q1, and its volumeis Vi1αβδ = Vt f i41 fi24 fi32. The vertices of the third tetrahed-ron are q1, qβ , qδ , qγ = f i31q3 + f i13q1, and its volume isVi1βδγ = Vt f i41 fi31 fi23. The locations of the vertices are depictedin figure 12(b). We obtainni = (Vi12αβ +V i1αβδ +Vi1βδγ)/Vt= f i42 fi32 + f i41 fi24 fi32 + f i41 fi31 fi23. (A17)Ji1 = [Vi12αβ +V i1αβδ(1+ f i14)+Vi1βδγ (1+ f i14 + f i13)]/(4niVt). (A18)Ji2 = [Vi12αβ (1+ f i24 + f i23)+V i1αβδ( fi24 + f i23)+Vi1βδγ fi23]/(4niVt). (A19)Ji3 = [Vi12αβ fi32 +V i1αβδ fi32 +Vi1βδγ ( fi32 + f i31)]/(4niVt).(A20)Ji4 = [Vi12αβ fi42 +V i1αβδ( fi42 + f i41)+Vi1βδγ fi41]/(4niVt).(A21)gi = 3( f i23 fi31 + f i32 fi24)/∆i41. (A22)I i1 = f i14/3+ f i13 fi31 fi23/(gi∆i41). (A23)I i2 = f i23/3+( f i24)2 f i32/(gi∆i41). (A24)I i3 = f i32/3+( f i31)2 f i23/(gi∆i41). (A25)I i4 = f i41/3+ f i42 fi24 fi32/(gi∆i41). (A26)A.4. ωi3 < ω < ωi4As shown in figure 12(c), the occupied part of the tet-rahedron is the full tetrahedron minus the tetrahedronwith vertices at q4, qβ = f i42q4 + f i24q2, qδ = f i41q4 + f i14q1and qγ = f i43q4 + f i34q3. Its volume is Vi4βδγ = Vt f i14 fi24 fi34.Therefore, the contribution of i-th tetrahedron to J(ω)is approximated as {Vt[F(q1)+F(q2)+F(q3)+F(q4)]−Vi4βδγ [F(q4)+F(qβ)+F(qγ)+F(qδ)]}/4. We obtainni = (1− f i14 fi24 fi34). (A27)Ji1 = [1− ( f i14)2 f i24 fi34]/(4ni). (A28)Ji2 = [1− f i14( fi24)2 f i34]/(4ni). (A29)Ji3 = [1− f i14 fi24( fi34)2]/(4ni). (A30)Ji4 = [1− f i14 fi24 fi34(1+ f i41 + f i42 + f i43)]/(4ni). (A31)gi = 3(1− ni)/(ωi4 −ω) = 3f i24 fi34/∆i41. (A32)I ik = f ik4/3, k= 1,2,3. (A33)I i4 = ( f i41 + f i42 + f i43)/3. (A34)A.5. ωi4 < ωThe tetrahedron is fully occupied, therefore,gi = 0. (A35)ni = 1. (A36)I ik = 0. (A37)J ik = 1/4. (A38)20J. Phys.: Condens. Matter 35 (2023) 353001 Topical ReviewAppendix B. Summation formulaeWhen the supercell lattice vectors are collinear to those of theprimitive cell,(as,bs,cs) = (ap,bp,cp)Mp→s. (B1)With Mp→s = diag(n1,n2,n3), it is proved in most textbooksthat the set of lattice vectorsRl = (ap,bp,cp)l1l2l3 , li = 0, . . . ,ni− 1 (B2)= (ap,bp,cp)l (B3)and their conjugated wave vectorsq= (a∗p ,b∗p ,c∗p)p1/n1p2/n2p3/n3 , pi = 0, . . . ,ni− 1 (B4)= (a∗s ,b∗s ,c∗s )p1p2p3 (B5)= (a∗s ,b∗s ,c∗s )p, (B6)fulfil the summation formulae1det(Mp→s)∑{Rl}eiq·Rl = δq,G, (B7)1det(Mp→s)∑{q}e−iq·Rl = δRl,RL , (B8)where G is any reciprocal lattice vector, and RL any supercelllattice vector.When the supercell is not collinear to the primitive cell, theset of lattice vectors, {Rl}, are those integer linear combin-ations of (ap,bp,cp) located within the supercell (as,bs,cs),and the set of wave vectors, {q}, are those integer linearcombinations of (a∗s ,b∗s ,c∗s ) located within the reciprocal cell(a∗p ,b∗p ,c∗p). Since the above summations can not easily be splitinto the product of three geometric series, it is not obviouswhether the above results are still valid.However, we have seen in section 3.1 that SNF like trans-formations can always be used to obtain primitive cell andsupercell which are collinear, and for which the summationformulae (B7) and (B8) are obviously valid. We have(ãp, b̃p, c̃p) = (ap,bp,cp)P−1, (B9)(ãs, b̃s, c̃s) = (as,bs,cs)Q, (B10)and therefore(ã∗p , b̃∗p , c̃∗p) = (a∗p ,b∗p ,c∗p)P⊺, (B11)(ã∗s , b̃∗s , c̃∗s ) = (a∗s ,b∗s ,c∗s )Q−⊺ (B12)= (ã∗p , b̃∗p , c̃∗p)D−1, (B13)while the set of lattice vectors and wave vectors,R̃l = (ãp, b̃p, c̃p)̃l, (B14)q̃= (ã∗s , b̃∗s , c̃∗s )p̃, (B15)with l̃i = 0, . . . ,Di− 1 and p̃i = 0, . . . ,Di− 1, fulfils,1det(D)∑{R̃l}ei q̃·R̃l = δq̃,G, (B16)1det(D)∑{q̃}e−i q̃·R̃l = δR̃l,RL. (B17)Since the same lattices are generated by the unimodulartransformations (B11) and (B12), we can writeR̃l = Rl+RL, (B18)q̃= q+G (B19)where RL is a supercell lattice vector used to bring R̃l withinthe original supercell (as,bs,cs), and G a reciprocal lat-tice vector used to bring q̃ to the original reciprocal cell(a∗p ,b∗p ,c∗p). Since q ·RL,G ·RL,G ·Rl are multiple of 2π, anddet(Mp→s) = det(D), we obtain equations (B7) and (B8) forarbitrary supercell and primitive cell related by an integer mat-rixMp→s. Notice however that the definition of the set of lat-tice vectors and wave vectors for which those summation for-mulae hold is now more general. The wave vectors are alongthe reciprocal vectors of the supercell and located within thereciprocal cell of the primitive lattice, while the lattice vectorsare along the primitive lattice vectors and located within thesupercell.Appendix C. Projection into primitive cell BlochstatesFrom the definition of the projection operator,T̂i P̂(q̃,q⋆)Ψlκ(q̃ν̃) = T̂i1NN∑j=1ei(q̃+q⋆)·Rlj T̂jΨlκ(q̃ν̃) (C1)= e−i(q̃+q⋆)·Rli1NN∑j=1ei(q̃+q⋆)·(Rlj+Rli )× T̂i T̂jΨlκ(q̃ν̃) (C2)21J. Phys.: Condens. Matter 35 (2023) 353001 Topical Review= e−i(q̃+q⋆)·Rli1NN∑j=1ei(q̃+q⋆)·(Rlj+Rli )× eiq̃·(Rlκ−Rlj−Rli )W(TiTj)−1(l)κ(q̃ν̃)(C3)= e−i(q̃+q⋆)·Rli1NN∑j=1eiq⋆·(Rlj+Rli )× eiq̃·RlκW(TiTj)−1(l)κ(q̃ν̃). (C4)Rlj +Rli may be outside the supercell, therefore wewriteRlj +Rli = Rlk +RL, whereRL is a supercell lattice vector. BecauseWlκ(q̃ν̃) has the periodicity of the supercell, and eiq⋆·RL = 1,we obtainT̂i P̂(q̃,q⋆)Ψlκ(q̃ν̃)= e−i(q̃+q⋆)·Rli1NN∑k=1eiq⋆·Rlk eiq̃·RlκWT−1k (l)κ(q̃ν̃) (C5)= e−i(q̃+q⋆)·Rli1NN∑k=1ei(q̃+q⋆)·Rlk eiq̃·(Rlκ−Rlk )WT−1k (l)κ(q̃ν̃)(C6)= e−i(q̃+q⋆)·Rli1NN∑k=1ei(q̃+q⋆)·Rlk T̂kΨlκ(q̃ν̃) (C7)= e−i(q̃+q⋆)·Rli P̂(q̃,q⋆)Ψlκ(q̃ν̃). (C8)ORCID iDsAtsushi Togo https://orcid.org/0000-0001-8393-9766Terumasa Tadano https://orcid.org/0000-0002-8132-2161References[1] Gonze X and Lee C 1997 Phys. Rev. B 55 10355[2] MacDonald A H, Vosko S H and Coleridge P T 1979 J. Phys.C: Solid State Phys. 12 2991[3] Blöchl P E, Jepsen O and Andersen O K 1994 Phys. Rev. B49 16223[4] Allen P B, Berlijn T, Casavant D A and Soler J M 2013 Phys.Rev. B 87 085322[5] Aroyo M I 2016 International Tables for Crystallography volA (Hoboken: Wiley)[6] Aroyo M I, Orobengoa D, de la Flor G, Tasci E S,Perez-Mato J M and Wondratschek H 2014 ActaCrystallogr. A 70 126[7] Togo A 2023 J. Phys. Soc. Japan 92 012001[8] Maradudin A A and Vosko S H 1968 Rev. Mod. Phys. 40 1[9] Shmueli U 2008 International Tables for Crystallography volB (Hoboken: Wiley)[10] Bradley T, Bradley C and Cracknell A 1972 The MathematicalTheory of Symmetry in Solids: Representation Theory forPoint Groups and Space Groups (Oxford: Clarendon)[11] Cracknell A P and Davies B L 1979 General Introduction andTables of Irreducible Representations of Space Groups(Kronecker Product Tables vol 1) (New York: IFI/Plenum)[12] Horton G K and Maradudin A A (eds) 1974 DynamicalProperties of Solids vol 1 (Amsterdam: North-HollandPublishing Company)[13] Niggli P 1928 Krystallographische und strukturtheoretischeGrundbegriffe Handbuch der Experimentalphysik vol7 (Leipzig: Akademische Verlagsgesellschaft)pp 108–76[14] Gruber B 1973 Acta Crystallogr. A 29 433[15] Křivý I and Gruber B 1976 Acta Crystallogr. A 32 297[16] Grosse-Kunstleve R W, Sauter N K and Adams P D 2004 ActaCrystallogr. A 60 1[17] Hart G L W and Forcade R W 2008 Phys. Rev. B 77 224115[18] Parlinski K, Li Z Q and Kawazoe Y 1997 Phys. Rev. Lett.78 4063[19] Wallace D C 1998 Thermodynamics of Crystals (New York:Dover Publications)[20] Togo A, Hayashi H, Tadano T, Tsutsui S and Tanaka I 2022 J.Phys.: Condens. Matter. 34 365401[21] Pick R M, Cohen M H and Martin R M 1970 Phys. Rev. B1 910[22] Giannozzi P, de Gironcoli S, Pavone P and Baroni S 1991Phys. Rev. B 43 7231[23] Gonze X, Charlier J-C, Allan D C and Teter M P 1994 Phys.Rev. B 50 13035[24] Hinuma Y, Pizzi G, Kumagai Y, Oba F and Tanaka I 2017Comput. Mater. Sci. 128 140[25] Togo A and Tanaka I 2018 Spglib: a software library forcrystal symmetry search (arXiv:1808.01590)[26] Monkhorst H J and Pack J D 1976 Phys. Rev. B 13 5188[27] Chadi D J and Cohen M L 1973 Phys. Rev. B 8 5747[28] Moreno J and Soler J M 1992 Phys. Rev. B 45 13891[29] Wisesa P, McGill K A and Mueller T 2016 Phys. Rev. B93 155109[30] Morgan W S, Jorgensen J J, Hess B C and Hart G L 2018Comput. Mater. Sci. 153 424[31] Morgan W S, Christensen J E, Hamilton P K, Jorgensen J J,Campbell B J, Hart G L and Forcade R W 2020 Comput.Mater. Sci. 173 109340[32] Chaput L, Togo A, Tanaka I and Hug G 2011 Phys. Rev. B84 094302[33] Togo A, Chaput L and Tanaka I 2015 Phys. Rev. B 91 094306[34] Alexei G H W, Maradudin A and Montroll E W (eds) 1963Theory of Lattice Dynamics in the HarmonicApproximation (Solid State Physics: Advances in Researchand Applications) (New York: Academic)[35] Monacelli L, Bianco R, Cherubini M, Calandra M, Errea Iand Mauri F 2021 J. Phys.: Condens. Matter33 363001[36] Errea I, Calandra M and Mauri F 2013 Phys. Rev. Lett.111 177002[37] Errea I, Calandra M and Mauri F 2014 Phys. Rev. B 89 064302[38] Bianco R, Errea I, Paulatto L, Calandra M and Mauri F 2017Phys. Rev. B 96 014111[39] Monacelli L, Errea I, Calandra M and Mauri F 2018 Phys. Rev.B 98 024106[40] Eriksson F, Fransson E and Erhart P 2019 Adv. Theory Simul.2 1800184[41] Tadano T and Tsuneyuki S 2015 Phys. Rev. B 92 054301[42] Wu Y-N, Saidi W A, Wuenschell J K, Tadano T, Ohodnicki P,Chorpening B and Duan Y 2020 J. Phys. Chem. Lett.11 2518[43] Tadano T and Tsuneyuki S 2018 J. Phys. Soc. Japan 87 041015[44] Moxon S et al 2022 J. Mater. Chem. A 10 186122https://orcid.org/0000-0001-8393-9766https://orcid.org/0000-0001-8393-9766https://orcid.org/0000-0002-8132-2161https://orcid.org/0000-0002-8132-2161https://doi.org/10.1103/PhysRevB.55.10355https://doi.org/10.1103/PhysRevB.55.10355https://doi.org/10.1088/0022-3719/12/15/008https://doi.org/10.1088/0022-3719/12/15/008https://doi.org/10.1103/PhysRevB.49.16223https://doi.org/10.1103/PhysRevB.49.16223https://doi.org/10.1103/PhysRevB.87.085322https://doi.org/10.1103/PhysRevB.87.085322https://doi.org/10.1107/S205327331303091Xhttps://doi.org/10.1107/S205327331303091Xhttps://doi.org/10.7566/JPSJ.92.012001https://doi.org/10.7566/JPSJ.92.012001https://doi.org/10.1103/RevModPhys.40.1https://doi.org/10.1103/RevModPhys.40.1https://doi.org/10.1107/S0567739473001063https://doi.org/10.1107/S0567739473001063https://doi.org/10.1107/S0567739476000636https://doi.org/10.1107/S0567739476000636https://doi.org/10.1107/S010876730302186Xhttps://doi.org/10.1107/S010876730302186Xhttps://doi.org/10.1103/PhysRevB.77.224115https://doi.org/10.1103/PhysRevB.77.224115https://doi.org/10.1103/PhysRevLett.78.4063https://doi.org/10.1103/PhysRevLett.78.4063https://doi.org/10.1088/1361-648X/ac7b01https://doi.org/10.1088/1361-648X/ac7b01https://doi.org/10.1103/PhysRevB.1.910https://doi.org/10.1103/PhysRevB.1.910https://doi.org/10.1103/PhysRevB.43.7231https://doi.org/10.1103/PhysRevB.43.7231https://doi.org/10.1103/PhysRevB.50.13035https://doi.org/10.1103/PhysRevB.50.13035https://doi.org/10.1016/j.commatsci.2016.10.015https://doi.org/10.1016/j.commatsci.2016.10.015https://arxiv.org/abs/1808.01590https://doi.org/10.1103/PhysRevB.13.5188https://doi.org/10.1103/PhysRevB.13.5188https://doi.org/10.1103/PhysRevB.8.5747https://doi.org/10.1103/PhysRevB.8.5747https://doi.org/10.1103/PhysRevB.45.13891https://doi.org/10.1103/PhysRevB.45.13891https://doi.org/10.1103/PhysRevB.93.155109https://doi.org/10.1103/PhysRevB.93.155109https://doi.org/10.1016/j.commatsci.2018.06.031https://doi.org/10.1016/j.commatsci.2018.06.031https://doi.org/10.1016/j.commatsci.2019.109340https://doi.org/10.1016/j.commatsci.2019.109340https://doi.org/10.1103/PhysRevB.84.094302https://doi.org/10.1103/PhysRevB.84.094302https://doi.org/10.1103/PhysRevB.91.094306https://doi.org/10.1103/PhysRevB.91.094306https://doi.org/10.1088/1361-648X/ac066bhttps://doi.org/10.1088/1361-648X/ac066bhttps://doi.org/10.1103/PhysRevLett.111.177002https://doi.org/10.1103/PhysRevLett.111.177002https://doi.org/10.1103/PhysRevB.89.064302https://doi.org/10.1103/PhysRevB.89.064302https://doi.org/10.1103/PhysRevB.96.014111https://doi.org/10.1103/PhysRevB.96.014111https://doi.org/10.1103/PhysRevB.98.024106https://doi.org/10.1103/PhysRevB.98.024106https://doi.org/10.1002/adts.201800184https://doi.org/10.1002/adts.201800184https://doi.org/10.1103/PhysRevB.92.054301https://doi.org/10.1103/PhysRevB.92.054301https://doi.org/10.1021/acs.jpclett.0c00183https://doi.org/10.1021/acs.jpclett.0c00183https://doi.org/10.7566/JPSJ.87.041015https://doi.org/10.7566/JPSJ.87.041015https://doi.org/10.1039/D1TA10072Fhttps://doi.org/10.1039/D1TA10072F Implementation strategies in phonopy and phono3py 1. Introduction 2. Crystal and symmetry 2.1. Crystal structure 2.2. Space group operation 2.3. Primitive cell 2.4. Reciprocal space 2.5. Phonon coordinates 2.6. BZ 3. Geometry of supercell model 3.1. Supercell construction 3.2. Commensurate points 4. Transformation between force constants and dynamical matrices 5. Non-analytical term correction 6. Regular grid in reciprocal primitive cell 6.1. Traditional regular grid 6.2. Generalized regular grid 6.3. Indexing of grid points 6.4. Symmetry of generalized regular grids 6.5. Double grid for subgrid shift 6.6. Grid points in BZ 6.7. Irreducible grid points 6.8. q point triplets 7. Tetrahedron method 7.1. Division into six tetrahedra 7.2. Functions 7.3. Integration weights 7.4. Use of symmetry of grid points 7.5. Triplets integration weights 8. Random displacement generation 9. Phonon band unfolding 10. Conclusion Appendix A. Functions used in the linear tetrahedron method A.1.  ωωi1 A.2.  ωi1 ωωi2 A.3.  ωi2 ωωi3 A.4.  ωi3 ωωi4 A.5.  ωi4 ω Appendix B. Summation formulae Appendix C. Projection into primitive cell Bloch states References