Group theoretical approach to computing phonons and their interactions

Phonons and their interactions are necessary for determining a wide range of materials properties. Here we present four independent advances which facilitate the computation of phonons and their interactions from first principles. First, we implement a group-theoretical approach to construct the order $\mathcal{N}$ Taylor series of a $d$-dimensional crystal purely in terms of space group irreducible derivatives (ID), which guarantees symmetry by construction and allows for a practical means of communicating and storing phonons and their interactions. Second, we prove that the smallest possible supercell which accommodates $\mathcal{N}$ given wave vectors in a $d$-dimensional crystal is determined using the Smith normal form of the matrix formed from the corresponding wave vectors; resulting in negligible computational cost to find said supercell, in addition to providing the maximum required multiplicity for uniform supercells at arbitrary $\mathcal{N}$ and $d$. Third, we develop a series of finite displacement methodologies to compute phonons and their interactions which exploit the first two developments: lone and bundled irreducible derivative (LID and BID) approaches. LID computes a single ID, or as few as possible, at a time in the smallest supercell possible, while BID exploits perturbative derivatives for some order less than $\mathcal{N}$ (e.g., Hellman-Feynman forces) in order to extract all ID in the smallest possible supercells using the fewest possible computations. Finally, we derive an equation for the order $N$ volume derivatives of the phonons in terms of the order $\mathcal{N}=N+2$ ID. Given that the former are easily computed, they can be used as a stringent, infinite ranged test of the ID. Our general framework is illustrated on graphene, yielding irreducible phonon interactions to fifth order. Additionally, we provide a cost analysis for the rocksalt structure at $\mathcal{N}=3$, demonstrating a massive speedup compared to popular finite displacement methods in the literature.

Figure 1

(a) A schematic of the structure of graphene. Yellow hexagons are lattice points, while circles represent carbon atoms. The first five smallest BvK supercells are pictured, and the lattice points of a corresponding FTG is given by Eq. (7). (b) The corresponding (color coded) five reciprocal lattice subcells which are repeated to tile the first Brillouin zone (FBZ); irreducible Brillouin zone is shaded grey. (c) Same as (b), but with the FBZ in the Wigner-Seitz cell convention. All $q$-points are labeled according to their star.

Figure 2

(a) Schematic of the graphene crystal structure, where yellow hexagons represent lattice points and circles represent carbon atoms; and each carbon atom is labeled by two integers which correspond to a translation in lattice coordinates. The ${\widehat{\mathbf{S}}}_K$ supercell is shown in red, and the corresponding WS cell is shown in blue and green for a centering on the first and second carbon atom, respectively. (b) Schematic showing $\mathbf{t}\in {\tilde{t}}_{\mathrm{BZ}}$ for ${\widehat{\mathbf{S}}}_{\mathrm{BZ}}={\widehat{\mathbf{S}}}_K$ along with the corresponding basis atoms. (c) Schematic showing how the basis atoms are translated back into the WS cell using some vector $\mathbf{t}{\widehat{\mathbf{S}}}_K\widehat{\mathbf{a}}$, with $\mathbf{t}\in {\mathbb{Z}}^d$, where the centering of the WS cell is on the first carbon atom. (d) Same as (c) but with the WS cell centered on the second carbon atom.

Figure 3

Phonons of graphene within DFT for ${\widehat{\mathbf{S}}}_{\mathrm{BZ}}=\widehat{\mathbf{1}}$, ${\widehat{\mathbf{S}}}_K$, and $12\widehat{\mathbf{1}}$, where data points are direct computational measurements and lines are Fourier interpolation. The irreducible derivatives for $\widehat{\mathbf{1}}$ and ${\widehat{\mathbf{S}}}_K$ are shown in Table 1; and the corresponding frequencies are obtained, in units of ${\mathrm{s}}^{-1}$, as ${\omega }_{\mathbf{q}}^{\alpha }=\sqrt{d_{\overline{\mathbf{q}}\mathbf{q}}^{\alpha \alpha }/m}$, where $m=12.011·1.0364\times {10}^{-28}\mathrm{eV}{\mathrm{s}}^2{\AA }^{-2}$ for carbon. The $y$ axis plots ${\omega }_{\mathbf{q}}^{\alpha }/\left(2\pi c\right)$, where $c$ is the speed of light in units of cm/s.

Figure 4

All panels display central finite difference calculations as a function of $\mathrm{\Delta }$; points are calculated values while the red line is a quadratic fit to a subset of points chosen by the algorithm defined in Sec. 3c. (a), (b), and (c) display a particular third-, fourth-, and fifth-order space group irreducible derivative, respectively, as obtained using the LID method with ${\mathrm{PD}}_1$. The horizontal dashed lines show the result of the SS-BID method, where many irreducible derivative are simultaneously extracted.

Figure 5

Complexity analysis for rocksalt structure at $\mathcal{N}=3$, comparing existing published methods (shengbte, phono3py, and aapl, taken from Ref. [85]) with our approaches (SS-BID and HS-BID); including (a) number of required DFT calculations and (b) time complexity assuming that the DFT calculations scale quadratically with system size [21, 22]. Existing methods only calculate the derivatives out to the $x$ neighbor shell (where $x$ is the horizontal axis) within supercell ${\widehat{\mathbf{S}}}_{\mathrm{BZ}}=x\widehat{\mathbf{1}}$, while our methods computes all derivatives within the corresponding supercell.

Figure 6

(a) Grüneisen parameters of graphene directly measured using identity strain derivatives of the phonons and Fourier interpolation (blue points and lines; diamonds and circles correspond to out-of-plane and in-plane modes, respectively) and via Eq. (50), which uses the cubic irreducible derivatives at various mesh densities. (b) and (c) follow the same conventions, but for the second and third strain derivatives, respectively; and both panels display the LID results for $2\widehat{\mathbf{1}}$ in gray, showing near perfect agreement with BID results.

Figure 7

A histogram of mean square error, associated with the quadratic fits to the central finite difference calculations as a function of $\mathrm{\Delta }$ within the SBB in the SS-BID approach, divided by the average magnitude of the SS-BID SBB derivatives. Note that values for $\mathcal{N}=2$ are multiplied by ten while the $\mathcal{N}=5$ values are divided by three for ease of viewing. FTG's of $6\widehat{\mathbf{1}}$, $2{\widehat{\mathbf{S}}}_K$, $2\widehat{\mathbf{1}}$, and $2\widehat{\mathbf{1}}$, are used for $\mathcal{N}=2\text{–}5$, respectively.

Figure 8

Comparison of the Grüneisen parameters in graphene obtained from cubic irreducible derivatives within ${\widehat{\mathbf{S}}}_{\mathrm{BZ}}=2{\widehat{\mathbf{S}}}_K$ using SS-BID. The blue curve uses our algorithm outlined in Sec. 3c to properly extrapolate $\mathrm{\Delta }$ to zero, while the other curves simply use a single value of $\mathrm{\Delta }$.

Figure 9

A plot comparing the irreducible derivatives computed using LID and SS-BID. The $x$ axis denotes relative magnitude $\left(\right|d_{\mathrm{LID}}|/\mathrm{RMS}\left(\left\{d_{\mathrm{LID}}\right\}\right))$ and the $y$ axis denotes relative error $\left(\right|d_{\mathrm{LID}}-d_{\mathrm{BID}}|/|d_{\mathrm{LID}}\left|\right)$. For $\mathcal{N}=3$ and 4, all irreducible derivatives are sampled, while only a subset is provided for $\mathcal{N}=5$.

Figure 10

Phonons of $\mathrm{ZrO}{}_2$ within DFT for ${\widehat{\mathbf{S}}}_{\mathrm{BZ}}=2{\widehat{\mathbf{S}}}_C$, where data points are direct computational measurements and lines are Fourier interpolation of the measurements. The irreducible derivatives for ${\widehat{\mathbf{S}}}_{\mathrm{BZ}}$ are shown in Table 3. When irreducible representations do not repeat at a given $\mathbf{q}$, the phonon frequency is given by ${\omega }_{\mathbf{q}}^{\alpha }=\sqrt{d_{\overline{\mathbf{q}}}^{\alpha }{}_{\mathbf{q}}^{\alpha }/m}$, where $m=m_i·1.0364\times {10}^{-28}\mathrm{eV}{\mathrm{s}}^2{\AA }^{-2}$, with either $m_O=15.9994$ or $m_{Zr}=91.224$. The $y$ axis plots ${\omega }_{\mathbf{q}}^{\alpha }\times {10}^{-12}/\left(2\pi \right)$, giving units of terahertz (THz). LO-TO splitting has not been incorporated.