 Harmonic Phonons Dynamical Matrix It is convenient to simplify the whole crystal to a supercell with three-dimensional boundary condition imposed on this supercell. This approach considerably reduces the number of explicitly treated atoms of a solid. The Phonon and PhononA softwares use the supercell approach. We remind that at the equilibrium all forces acting on supercell atoms vanish. The simplified version of the computed method relies on determinding forces F appearing on all supercell atoms due to a displacement U of a single atom, one at the time, and calculating from them the crystal force constants Φ. The number of atomic displacements can be highly reduced when taking into account the crystal symmetry.  In the nearby figures for the cubic NiAl intermetalic only two displacements along single direction Z are needed, one for Al and second for Ni atom. All other displacements are symmetrically equivalent to these two mentioned. The forces F are called Hellmann-Feynman forces. The potential energy of the crystal can be expanded in terms of atomic displacements from equilibrium positions. In harmonic approximation the expanssion is truncated at the second-order term, and then it looks like: V = V(.... U(n, μ),...)  =  V(....0,...)   +  (½)Σn,μ,m,ν    Φ(n,μ,m,ν)  U(n,μ) U(m,ν)  +  ...      (1) V = V(.... U(n, μ),...)  +  (½) Σn,μ,m,ν    Φ(n,μ,m,ν)  U(n,μ) U(m,ν)  +  ...      (1) The force acting on an atom (n, μ) is: F(n,μ) = - ∂V(.... U(n, μ),...) / ∂ U(n,ν)       (2) Thus, using Eq.(1) we find relation between forces and force constants as: F(n,μ) = - Σm,ν Φ(n,μ,m,ν)  U(m,ν)       (3) Hence, the force acting on an atom (n, μ) arises from displacements U(m,ν) of all supercell atoms (m,ν). Eq.(3) is the master equation of the direct method. And for simplicity one applies only a single atomic displacement. This displacement U(m,ν) is an input for the ab initio calculations, while the Hellmann-Feynman forces F(n,μ) are found as a result of a single ionic loop of optimization in ab initio calculations. The master equation seems to be a linear algebraic equation. In fact, the solution of the master equation (3) is far from being simple, since the number of known and unknown variables is different, and the space group of the crystal requires to impose on Eq. (3) additional symmetry constrains . While finding the force constants, the dynamical matrix is straightforward to construct. Its elements are determined as: D(k;0,μ,m,ν) = 1 /(√M0,μ √ Mm,ν) Σm,ν Φ(0,μ,m,ν)  exp[ -ik(R(0,μ) - R(m,ν)]   (4) The phonon frequencies and polarization vectors are now obtained from the solution of the eigenvalue problem ω2(k,j)  E(k,j) =  D(k)  E(k,j)       (5) To obtain the reliable phonon dispersion relations the supercell size should be about 6 - 8 Å. The magnitude of the force constants beyond that distance is expected to be negligible. In such a case the phonon dispersion relations can be calculated exactly. But even if the supercell size is small, exact phonon frequencies can be obtained for the wave vectors commensurate with the supercell size. Then, the phonon dispersion curves are the symmetry controlled interpolations between the exact points. To write a software which handles the Eqs.(3-5) for all possible crystals and structure needed a consistent and convenient formulation of these equations. An algoritm in matrix formulation was proposed. It can find force constants from interatomic forces, so called Hellmann-Feynman forces, and in addition it can take into account the symmetry elements of the crystal space groups, and the necessary invariants. It has been published in the Letter of K.Parlinski, Z.Q.Li and Y.Kawazoe , and written with all the details by Parlinski in the Manual of the first addition of Phonon ver. 1.01 in 1998 . Reference  K.Parlinski, Z.Q.Li and Y.Kawazoe, First-principle determination of the soft mode in cubic ZrO2, Phys.Rev.Lett. 78, 4063 (1997).  K.Parlinski, Phonon Manual, ver. 1.04, Cracow (1998).