# Projector-based efficient estimation of force constants

Atsuto Seko<EMAIL_ADDRESS>Department of Materials Science and
Engineering, Kyoto University, Kyoto 606-8501, Japan  Atsushi Togo  Center for
Basic Research on Materials National Institute for Materials Science, Tsukuba,
Ibaraki 305-0047, Japan

###### Abstract

Estimating force constants for crystal structures is a fundamental aspect of
calculating various phonon-related properties. However, it can be a
challenging task, particularly when dealing with a large number of atoms or
when third and higher-order force constants are needed. In this study, we
propose an efficient approach that involves constructing a complete
orthonormal basis set of the force constants. We formulate this approach using
projector matrices and their eigenvectors corresponding to the requirements
for the force constants. This basis set allows us to infer the force constants
from displacement-force datasets precisely. Our efficient algorithms for the
basis-set construction and force constant estimation from a displacement-force
dataset are implemented in the symfc code Seko and Togo . Furthermore, several
applications are demonstrated in this study, indicating that the current
approach helps to obtain force constants efficiently and precisely.

##  I Introduction

Force constants of crystal structures play a crucial role in calculating
various phonon-related properties Dove (1993); Wallace (1998); Ziman (1960).
To estimate these force constants, we typically use supercell approaches in
conjunction with the finite displacement method Parlinski  _et al._ (1997);
Hellman  _et al._ (2013); Tadano and Tsuneyuki (2018); Eriksson  _et al._
(2019); Togo (2023). Moreover, in self-consistent phonon approaches, effective
force constants are often determined from supercell structures that replicate
the atomic distribution at a specific temperature Errea  _et al._ (2013); van
Roekeghem  _et al._ (2021). In these approaches, random or systematic atomic
displacements are introduced into supercells. The forces acting on atoms in
the supercells are calculated, and then supercell force constants are
estimated from the dataset composed of the displacements and forces. However,
determining the force-constant elements necessitates numerous supercells with
different configurations of displacements because the force-constant elements
that need to be determined are plentiful.

The force-constant elements that need to be determined can be reduced into a
basis set for representing the force constants by applying permutation
symmetry rules, which are derived from the definition of the force constants
as the derivatives of the potential energy, and sum rules, which are derived
from the invariance with infinitesimal translations Wallace (1998). In
addition, symmetric properties of the structure can also reduce the force-
constant elements Tadano and Tsuneyuki (2018); Eriksson  _et al._ (2019); Nye
(1985). A significantly reduced but complete orthonormal basis set allows for
efficient and precise calculation of force constants that meet all the
requirements for the force constants El-Batanouny and Wooten (2008);
Venkataraman  _et al._ (1975). However, constructing the basis set and
estimating the force constants using the basis set can be computationally
demanding, especially when dealing with many atoms or when third and higher-
order force constants are needed.

Various approaches have been developed to estimate force constants precisely
using compressed sensing, as reported in the literature Zhou  _et al._ (2014);
Tadano and Tsuneyuki (2015). Although there are other approaches available
Tadano and Tsuneyuki (2018); Eriksson  _et al._ (2019); Hellman  _et al._
(2013), the force constant estimation remains computationally challenging. In
this study, we propose a projector-based efficient approach that involves
constructing a complete orthonormal basis set of force constants. We introduce
several compression techniques that allow the basis set to be efficiently
reduced. These techniques compress a projection matrix using eigenvectors of
another projection matrix onto subspaces defined by partial requirements for
the force constants. We have implemented the current procedures for efficient
basis-set construction and force constant estimation from a displacement-force
dataset in the symfc code Seko and Togo . We also demonstrate various
applications of the current approach, indicating that it leads to obtaining
force constants efficiently and precisely.

##  II Formulation of projection matrices

###  II.1 Crystal potential and force constants

Force constants of a crystal are derived from the Taylor expansion of the
potential energy $\mathcal{V}$ with respect to atomic displacements. In this
study, we consider supercell force constants, which are estimated from a set
of finite displacements and forces acting on atoms in supercells. For a
supercell with $N$ atoms, the potential energy is represented by

$\displaystyle\mathcal{V}=\Theta_{0}$ $\displaystyle+$
$\displaystyle\sum_{i,\alpha}\Theta_{i\alpha}u_{i\alpha}$ $\displaystyle+$
$\displaystyle\frac{1}{2}\sum_{i\alpha,j\beta}\Theta_{i\alpha,j\beta}u_{i\alpha}u_{j\beta}$
$\displaystyle+$
$\displaystyle\frac{1}{3!}\sum_{i\alpha,j\beta,k\gamma}\Theta_{i\alpha,j\beta,k\gamma}u_{i\alpha}u_{j\beta}u_{k\gamma}$
$\displaystyle+$ $\displaystyle\cdots,$

where $i$, $j$, and $k$ are the atomic indices within the supercell, and
$\alpha$, $\beta$, and $\gamma$ are the Cartesian components. The second-order
and third-order supercell force constants are represented by
$\Theta_{i\alpha,j\beta}$ and $\Theta_{i\alpha,j\beta,k\gamma}$, respectively.
In addition, we assume that the first-order force constants $\Theta_{i\alpha}$
is equal to zero because atomic displacements are measured from their
equilibrium positions.

###  II.2 Projection matrix for second-order force constants

In the current formulation, the second-order force constants are represented
in a column vector form with $9N^{2}$ elements as

$\bm{\theta}^{\rm(FC2)}=[\Theta_{1x,1x},\Theta_{1x,1y},\Theta_{1x,1z},\Theta_{1x,2x},\cdots]^{\top},$
(2)

which are listed in the ascending order of $(i,\alpha,j,\beta)$. We start by
considering force constants that are invariant with space group operations
$\mathcal{G}$ for the target structure. As we are discussing supercell force
constants, the lattice translation subgroup in $\mathcal{G}$ is limited to
those of the lattice points inside the supercell. An operation
$\hat{g}\in\mathcal{G}$ acts on the force constant vector as

$\hat{g}\bm{\theta}^{\rm(FC2)}=\bm{\theta}^{\rm(FC2)}\left[\bm{\Gamma}(\hat{g})\otimes\bm{\Gamma}(\hat{g})\right],$
(3)

where $\bm{\Gamma}(\hat{g})$ denotes the matrix representation of operation
$\hat{g}$ for the first-order force constants $\\{\Theta_{i\alpha}\\}$. We use
the group-theoretical projection operator method El-Batanouny and Wooten
(2008) to reduce the Kronecker product of representations to the identity
irreducible representation. The projection matrix for the identity irreducible
representation is given as

$\bm{P}_{\rm
spg}=\frac{1}{|\mathcal{G}|}\sum_{\hat{g}\in\mathcal{G}}\bm{\Gamma}(\hat{g})\otimes\bm{\Gamma}(\hat{g}),$
(4)

where $|\mathcal{G}|$ denotes the order of group $\mathcal{G}$.

Then, we introduce projection matrices onto force constants that satisfy the
sum rules and permutation symmetry rules. The sum rule for $(i\alpha,\beta)$
and the permutation symmetry rule for $(i\alpha,j\beta)$ are expressed as

$\sum_{j}\Theta_{i\alpha,j\beta}=0$ (5)

and

$\Theta_{i\alpha,j\beta}-\Theta_{j\beta,i\alpha}=0,$ (6)

respectively. These rules are independent linear constraints for the force
constants. We can represent them in a matrix form as
$\bm{C}^{\top}\bm{\theta}^{\rm(FC2)}=0$. Therefore, a force-constant basis set
that satisfies these linear constraints spans the complementary subspace of
the kernel of $\bm{C}^{\top}$ Yanai  _et al._ (2011). The kernel of
$\bm{C}^{\top}$ is written as

${\rm
Ker}(\bm{C}^{\top})=\\{\bm{\theta}^{\rm(FC2)}\in\mathbb{R}^{9N^{2}}|\bm{C}^{\top}\bm{\theta}^{\rm(FC2)}=0\\},$
(7)

and the orthogonal projection matrix onto the complementary subspace of ${\rm
Ker}(\bm{C}^{\top})$ is given by Strang (2023); Yanai  _et al._ (2011)

$\bm{P}_{\rm perm\cap sum}=\bm{I}-\bm{C}(\bm{C^{\top}C})^{-1}\bm{C}^{\top},$
(8)

where $\bm{I}$ denotes the unit matrix of size $9N^{2}$.

The column unit vector $\bm{v}_{i\alpha\beta}$ represents the sum rule for
$(i,\alpha,\beta)$, which can be expressed as
$\bm{v}_{i\alpha\beta}^{\top}\bm{\theta}^{\rm(FC2)}=0$, in accordance with Eq.
(5). This means that the vector $\bm{v}_{i\alpha\beta}$ has only $N$ non-zero
elements. Moreover, the set of vectors $\\{\bm{v}_{i\alpha\beta}\\}$ are
mutually orthonormal. Similarly, the column unit vector
$\bm{w}_{i\alpha,j\beta}$ describes the permutation symmetry rule for the
combination of $i\alpha$ and $j\beta$ as
$\bm{w}_{i\alpha,j\beta}^{\top}\bm{\theta}^{\rm(FC2)}=0$. Therefore, the
vector $\bm{w}_{i\alpha,j\beta}$ has only two non-zero elements, which are
$1/\sqrt{2}$ and $-1/\sqrt{2}$. Additionally, the set of vectors
$\\{\bm{w}_{i\alpha,j\beta}\\}$, with $i\alpha<j\beta$ in the lexicographic
order, is orthonormal.

Using the complete set of $\bm{v}_{i\alpha\beta}$ and
$\bm{w}_{i\alpha,j\beta}$, matrix $\bm{C}$ is described as

$\bm{C}=[\bm{C}_{v},\bm{C}_{w}],$ (9)

where

$\bm{C}_{v}=[\cdots,\bm{v}_{i\alpha\beta},\cdots],\>\>\bm{C}_{w}=[\cdots,\bm{w}_{i\alpha,j\beta},\cdots].$
(10)

Note that two vectors $\bm{v}_{i\alpha\beta}$ and $\bm{w}_{i\alpha,j\beta}$
are linearly independent, without necessarily being orthogonal. Therefore, the
projection matrix onto the complementary subspace of ${\rm
Ker}(\bm{C}^{\top})$ requires an inverse matrix computation as given by Eq.
(8).

Finally, we present the full projection matrix onto the intersection of vector
subspaces determined by projection matrices $\bm{P}_{\rm spg}$ and
$\bm{P}_{\rm perm\cap sum}$. Since projection matrices $\bm{P}_{\rm spg}$ and
$\bm{P}_{\rm perm\cap sum}$ are commutative as shown in the Appendix B, the
full projection matrix can be written as $\bm{P}=\bm{P}_{\rm spg}\bm{P}_{\rm
perm\cap sum}$ Yanai  _et al._ (2011); Piziak  _et al._ (1999). By solving the
eigenvalue problem for the full projection matrix, i.e.,

$\bm{P}\bm{b}=\bm{b},$ (11)

we can obtain an orthonormal basis set for representing force constants. This
basis set is invariant with space group operations and satisfies the sum rules
and permutation symmetry rules. Since all non-zero eigenvalues of the
projection matrix are one, the number of independent basis vectors is equal to
the trace of the projection matrix.

###  II.3 Projection matrix for third-order force constants

Third-order force constants are represented in a vector form with $27N^{3}$
elements as

$\bm{\theta}^{\rm(FC3)}=[\Theta_{1x,1x,1x},\Theta_{1x,1x,1y},\Theta_{1x,1x,1z},\Theta_{1x,1x,2x},\cdots]^{\top},$
(12)

which are listed in the ascending order of $(i,\alpha,j,\beta,k,\gamma)$. An
operation $\hat{g}\in\mathcal{G}$ acts on the third-order force constant
vector as

$\hat{g}\bm{\theta}^{\rm(FC3)}=\bm{\theta}^{\rm(FC3)}\left[\bm{\Gamma}(\hat{g})\otimes\bm{\Gamma}(\hat{g})\otimes\bm{\Gamma}(\hat{g})\right].$
(13)

Therefore, the projection matrix for the identity irreducible representation
is given as

$\bm{P}_{\rm
spg}=\frac{1}{|\mathcal{G}|}\sum_{\hat{g}\in\mathcal{G}}\bm{\Gamma}(\hat{g})\otimes\bm{\Gamma}(\hat{g})\otimes\bm{\Gamma}(\hat{g}).$
(14)

The sum rules and permutation symmetry rules can also be represented in a
matrix form as $\bm{C}^{\top}\bm{\theta}^{(\rm FC3)}=0$, which is the same as
the representation of the rules used for second-order force constants. The
kernel of $\bm{C}^{\top}$ is written as

${\rm Ker}(\bm{C}^{\top})=\\{\bm{\theta}^{(\rm
FC3)}\in\mathbb{R}^{27N^{3}}|\bm{C}^{\top}\bm{\theta}^{(\rm FC3)}=0\\},$ (15)

and the orthogonal projection matrix onto the complementary subspace of ${\rm
Ker}(\bm{C}^{\top})$ is given by Strang (2023); Yanai  _et al._ (2011)

$\bm{P}_{\rm perm\cap sum}=\bm{I}-\bm{C}(\bm{C^{\top}C})^{-1}\bm{C}^{\top},$
(16)

where $\bm{I}$ denotes the unit matrix of size $27N^{3}$.

The sum rule for $(i,\alpha,j,\beta,\gamma)$ is written as

$\sum_{k}\Theta_{i\alpha,j\beta,k\gamma}=0,$ (17)

hence, column unit vector $\bm{v}_{i{\alpha}j\beta\gamma}$ represents the sum
rule for $(i,\alpha,j,\beta,\gamma)$, which can be expressed as
$\bm{v}_{i{\alpha}j\beta\gamma}^{\top}\bm{\theta}^{(\rm FC3)}=0$. The vector
$\bm{v}_{i{\alpha}j\beta\gamma}$ has only $N$ non-zero elements among
$27N^{3}$ elements, and the set of vectors
$\\{\bm{v}_{i{\alpha}j\beta\gamma}\\}$ are mutually orthonormal.

Representing the permutation symmetry rules is more complicated than applied
to the second-order force constants. The permutation symmetry rule for
$(i\alpha,j\beta,k\gamma)$ is expressed as

$\displaystyle\Theta_{i\alpha,j\beta,k\gamma}-\Theta_{i\alpha,k\gamma,j\beta}$
$\displaystyle=$ $\displaystyle 0$
$\displaystyle\Theta_{i\alpha,j\beta,k\gamma}-\Theta_{j\beta,i\alpha,k\gamma}$
$\displaystyle=$ $\displaystyle 0$
$\displaystyle\Theta_{i\alpha,j\beta,k\gamma}-\Theta_{j\beta,k\gamma,i\alpha}$
$\displaystyle=$ $\displaystyle 0$ (18)
$\displaystyle\Theta_{i\alpha,j\beta,k\gamma}-\Theta_{k\gamma,i\alpha,j\beta}$
$\displaystyle=$ $\displaystyle 0$
$\displaystyle\Theta_{i\alpha,j\beta,k\gamma}-\Theta_{k\gamma,j\beta,i\alpha}$
$\displaystyle=$ $\displaystyle 0,$

when all three indices $i\alpha$, $j\beta$, and $k\gamma$ are different. These
equations can be described by a matrix composed of the column unit vectors
$\bm{W}_{i\alpha,j\beta,k\gamma}$ as
$\bm{W}_{i\alpha,j\beta,k\gamma}^{\top}\bm{\theta}^{(\rm FC3)}=0$. Each column
vector of the matrix $\bm{W}_{i\alpha,j\beta,k\gamma}$ has two non-zero
elements, which are $1/\sqrt{2}$ and $-1/\sqrt{2}$. The column vectors are not
orthogonal but linearly independent. The permutation symmetry rule for
$(i\alpha,i\alpha,k\gamma)$ is given as

$\displaystyle\Theta_{i\alpha,i\alpha,k\gamma}-\Theta_{i\alpha,k\gamma,i\alpha}$
$\displaystyle=$ $\displaystyle 0$
$\displaystyle\Theta_{i\alpha,i\alpha,k\gamma}-\Theta_{k\gamma,i\alpha,i\alpha}$
$\displaystyle=$ $\displaystyle 0.$ (19)

In this case, the matrix composed of column unit vectors
$\bm{W}_{i\alpha,i\alpha,k\gamma}$ is also written as
$\bm{W}_{i\alpha,i\alpha,k\gamma}^{\top}\bm{\theta}^{(\rm FC3)}=0$, where each
column vector has two non-zero elements, which are $1/\sqrt{2}$ and
$-1/\sqrt{2}$.

Using a linearly-independent complete set of $\bm{v}_{i{\alpha}j\beta\gamma}$
and $\bm{W}_{i\alpha,j\beta,k\gamma}$, matrix $\bm{C}$ is described as

$\bm{C}=[\bm{C}_{v},\bm{C}_{w}],$ (20)

where

$\bm{C}_{v}=[\cdots,\bm{v}_{i{\alpha}j\beta\gamma},\cdots],\>\>\bm{C}_{w}=[\cdots,\bm{W}_{i\alpha,j\beta,k\gamma},\cdots].$
(21)

Finally, the full projection matrix is given as $\bm{P}=\bm{P}_{\rm
spg}\bm{P}_{\rm perm\cap sum}$ Yanai  _et al._ (2011); Piziak  _et al._
(1999). By solving the eigenvalue problem for the full projection matrix,
$\bm{P}\bm{b}=\bm{b}$, an orthonormal basis set for representing third-order
force constants is obtained.

##  III Efficient basis set construction

###  III.1 Decomposition of projection matrix

To efficiently obtain the projection matrix of $\bm{P}_{\rm perm\cap sum}$ and
formulate a more efficient procedure, we use an approximate projection matrix
instead of computing the inverse of $\bm{C}^{\top}\bm{C}$, which can be
computationally prohibitive. This involves considering individual projection
matrices for the sum rules and permutation symmetry rules. For the second-
order force constants, orthogonal projection matrices $\bm{P}_{\bm{v}}$ and
$\bm{P}_{\bm{w}}$ onto the complementary subspaces of ${\rm
Ker}(\bm{C}_{v}^{\top})$ and ${\rm Ker}(\bm{C}_{w}^{\top})$ are expressed
using orthonormal vector sets $\\{\bm{v}_{i\alpha\beta}\\}$ and
$\\{\bm{w}_{i\alpha,j\beta}\\}$ as Strang (2023)

$\displaystyle\bm{P}_{v}$ $\displaystyle=$
$\displaystyle\bm{I}-\sum_{i=1}^{N}\sum_{\alpha=1}^{3}\sum_{\beta=1}^{3}\bm{v}_{i\alpha\beta}\bm{v}_{i\alpha\beta}^{\top}$
(22) $\displaystyle=$ $\displaystyle\bm{I}-\bm{C}_{v}\bm{C}_{v}^{\top}$

and

$\displaystyle\bm{P}_{w}$ $\displaystyle=$
$\displaystyle\bm{I}-\sum_{(i\alpha,j\beta)}\bm{w}_{i\alpha,j\beta}\bm{w}_{i\alpha,j\beta}^{\top}$
(23) $\displaystyle=$ $\displaystyle\bm{I}-\bm{C}_{w}\bm{C}_{w}^{\top},$

respectively. For the third-order force constants, orthogonal projection
matrices $\bm{P}_{\bm{v}}$ and $\bm{P}_{\bm{w}}$ onto the complementary
subspaces of ${\rm Ker}(\bm{C}_{v}^{\top})$ and ${\rm Ker}(\bm{C}_{w}^{\top})$
are given as

$\displaystyle\bm{P}_{v}$ $\displaystyle=$
$\displaystyle\bm{I}-\sum_{i=1}^{N}\sum_{\alpha=1}^{3}\sum_{j=1}^{N}\sum_{\beta=1}^{3}\sum_{\gamma=1}^{3}\bm{v}_{i\alpha
j\beta\gamma}\bm{v}_{i\alpha j\beta\gamma}^{\top}$ (24) $\displaystyle=$
$\displaystyle\bm{I}-\bm{C}_{v}\bm{C}_{v}^{\top}$

and

$\displaystyle\bm{P}_{w}$ $\displaystyle=$
$\displaystyle\bm{I}-\bm{C}_{w}(\bm{C}_{w}^{\top}\bm{C}_{w})^{-1}\bm{C}_{w}^{\top},$
(25)

respectively.

Using the non-commutative matrices $\bm{P}_{v}$ and $\bm{P}_{w}$, the
projection matrix $\bm{P}_{\rm perm\cap sum}$ is expressed using the von
Neumann–Halperin theorem as von Neumann (1949); Halperin (1962); Piziak  _et
al._ (1999); Galántai (2013)

$\bm{P}_{\rm perm\cap
sum}=\lim_{p\to\infty}\bm{P}_{w}(\bm{P}_{v}\bm{P}_{w})^{p}.$ (26)

To obtain an approximation of the projection matrix $\bm{P}_{\rm perm\cap
sum}$, a finite value of $p$ can be used. However, the matrix product
$\bm{P}_{\rm perm\cap sum}\simeq\bm{P}_{w}(\bm{P}_{v}\bm{P}_{w})^{1}$, i.e.,
$p=1$, suffices to derive a basis set for the full projection matrix. By
solving the eigenvalue problem of the approximated full projection matrix
given as

$\tilde{\bm{P}}=\bm{P}_{\rm spg}\bm{P}_{w}\bm{P}_{v}\bm{P}_{w}$ (27)

and removing its eigenvectors with eigenvalues less than one, we can obtain
the basis set for the full projection matrix.

###  III.2 Matrix compression

When computing the product of two large projection matrices of size
$(9N^{2},9N^{2})$ for the second order or $(27N^{3},27N^{3})$ for the third
order, it can be computationally expensive. To avoid allocating large matrices
and computing the products of large matrices, we use a compression technique
that involves the eigenvalue decomposition of projection matrices. In the
eigenvalue decomposition, a projection matrix $\bm{P}_{i}$ is decomposed by
its eigenvectors $\bm{B}_{i}$ as

$\bm{P}_{i}=\bm{B}_{i}\bm{B}_{i}^{\top},$ (28)

where $\bm{B}_{i}$ is of size $(9N^{2},n_{i})$ or $(27N^{3},n_{i})$. If the
eigenvectors of the projection matrix are well-compressed, sparse, and
efficiently obtained, they can be helpful in compressing large matrices.

Let us consider the case where product $\bm{P}_{j}\bm{P}_{i}\bm{P}_{j}$ can be
included in a product representation of the full projection matrix. If
$\bm{P}_{i}$ and $\bm{P}_{j}$ are commutative, i.e.,
$[\bm{P}_{i},\bm{P}_{j}]=0$, product $\bm{P}_{j}\bm{P}_{i}\bm{P}_{j}$ is equal
to $\bm{P}_{i}\bm{P}_{j}$ as

$\bm{P}_{j}\bm{P}_{i}\bm{P}_{j}=\bm{P}_{i}\bm{P}_{j}^{2}=\bm{P}_{i}\bm{P}_{j},$
(29)

which is derived using the idempotent property of the projection matrix,
$\bm{P}_{j}^{2}=\bm{P}_{j}$. When eigenvectors of $\bm{P}_{j}$, $\bm{B}_{j}$,
are known and well compressed, $\bm{P}_{i}$ may be efficiently compressed by
$\bm{B}_{j}$, as

$\displaystyle\bm{P}_{j}\bm{P}_{i}\bm{P}_{j}$ $\displaystyle=$
$\displaystyle\bm{B}_{j}\bm{B}_{j}^{\top}\bm{P}_{i}\bm{B}_{j}\bm{B}_{j}^{\top}\>\>\>(\bm{P}_{j}=\bm{B}_{j}\bm{B}_{j}^{\top})$
(30) $\displaystyle=$ $\displaystyle\bm{B}_{j}\bm{P}_{i\downarrow
j}\bm{B}_{j}^{\top},$

where $\bm{P}_{i\downarrow j}$ denotes the compressed matrix of $\bm{P}_{i}$
by basis set $\bm{B}_{j}$ as

$\bm{P}_{i\downarrow j}=\bm{B}_{j}^{\top}\bm{P}_{i}\bm{B}_{j}.$ (31)

Using this compression, the eigenvectors of $\bm{P}_{j}\bm{P}_{i}\bm{P}_{j}$
can be efficiently calculated as follows. (1) If eigenvectors of
$\bm{P}_{i}=\bm{B}_{i}\bm{B}_{i}^{\top}$ are derived without constructing
$\bm{P}_{i}$ and solving the eigenvalue problem for $\bm{P}_{i}$ as shown
below, $\bm{P}_{i\downarrow j}$ is written as

$\bm{P}_{i\downarrow
j}=\bm{B}_{j}^{\top}\bm{B}_{i}\bm{B}_{i}^{\top}\bm{B}_{j}.$ (32)

In general, the size of $\bm{B}_{j}^{\top}\bm{B}_{i}$ is small, so we firstly
compute $\bm{A}=\bm{B}_{j}^{\top}\bm{B}_{i}$, and then compute
$\bm{P}_{i\downarrow j}=\bm{A}\bm{A}^{\top}$. If it is necessary to construct
$\bm{P}_{i}$, we compute $\bm{P}_{i\downarrow j}$ using Eq. (31). (2) The
eigenvalue problem of compressed projection matrix $\bm{P}_{i\downarrow j}$ is
solved as $\bm{P}_{i\downarrow j}=\bm{B}_{i\downarrow j}\bm{B}_{i\downarrow
j}^{\top}$, where $\bm{B}_{i\downarrow j}$ is the compressed eigenvectors by
$\bm{B}_{j}$. (3) The eigenvectors of $\bm{P}_{j}\bm{P}_{i}\bm{P}_{j}$,
$\bm{B}_{i{\cdot}j}$, is derived as

$\bm{B}_{i{\cdot}j}=\bm{B}_{j}\bm{B}_{i\downarrow j},$ (33)

and $\bm{P}_{j}\bm{P}_{i}\bm{P}_{j}$ is rewritten by

$\displaystyle\bm{P}_{j}\bm{P}_{i}\bm{P}_{j}$ $\displaystyle=$
$\displaystyle\bm{B}_{i{\cdot}j}\bm{B}_{i{\cdot}j}^{\top}$ (34)
$\displaystyle=$ $\displaystyle\bm{B}_{j}\bm{B}_{i\downarrow
j}\bm{B}_{i\downarrow j}^{\top}\bm{B}_{j}^{\top}.$

This matrix compression procedure is frequently used in the reduction of the
full projection matrix as shown below.

###  III.3 Application of sum rules

We utilize the equation of $\bm{P}_{v}=\bm{I}-\bm{C}_{v}\bm{C}_{v}^{\top}$,
which is presented in Eqs. (22) and (24), to reduce the full projection matrix
by the sum rules. This method is efficient because both the unit matrix and
$\bm{C}_{v}$ are sparse matrices. On the other hand, the eigenvectors of
$\bm{P}_{v}$, $\bm{B}_{v}$, typically contains a larger number of non-zero
elements than $\bm{C}_{v}$.

###  III.4 Application of permutation symmetry rules

The eigenvectors of the projection matrix onto the complementary space of
${\rm Ker}(\bm{C}_{w}^{\top})$, $\bm{P}_{w}$, can be efficiently obtained
without actually constructing the projection matrix $\bm{P}_{w}$. These
eigenvectors form a basis for the vector space that satisfies the permutation
symmetry rules. For the second order, the permutation symmetry rule for
$(i\alpha,j\beta)$ is given by Eq. (23), and is represented by a single vector
$\bm{w}_{i\alpha,j\beta}$ as $\bm{w}_{i\alpha,j\beta}^{\top}\bm{\theta}=0$.
Because the set of vectors $\\{\bm{w}_{i\alpha,j\beta}\\}$ is orthonormal, the
vector $\bm{w}_{i\alpha,j\beta}^{\perp}$ orthogonal to
$\bm{w}_{i\alpha,j\beta}$ can be an eigenvector of $\bm{P}_{w}$, which has
only two non-zero elements of $1/\sqrt{2}$ for columns $(i,\alpha,j,\beta)$
and $(j,\beta,i,\alpha)$. If $i\alpha=j\beta$, no permutation symmetry rules
should be considered. Therefore, $\bm{w}_{i\alpha,i\alpha}^{\perp}$ can be
defined as the vector that has only a non-zero element of one for column
$(i,\alpha,i,\alpha)$. Thus, the eigenvectors of $\bm{P}_{w}$, $\bm{B}_{w}$,
is written as

$\bm{B}_{w}=[\cdots,\bm{w}_{i\alpha,j\beta}^{\perp},\cdots]^{\top},$ (35)

where vectors $\\{\bm{w}_{i\alpha,j\beta}^{\perp}\\}$ are also mutually
orthonormal. The size of $\bm{B}_{w}$ is $(9N^{2},n_{w})$, where $n_{w}\simeq
9N^{2}/2$. The number of non-zero elements is $9N^{2}$ in $\bm{B}_{w}$, which
is much smaller than the total number of elements in $\bm{B}_{w}$.

For the third order, the permutation symmetry rule for
$(i\alpha,j\beta,k\gamma)$ can be found in Eq. (II.3), which consists of five
equations. These equations are represented in a matrix form as
$\bm{W}_{i\alpha,j\beta,k\gamma}^{\top}\bm{\theta}=0$, and the submatrix of
$\bm{W}_{i\alpha,j\beta,k\gamma}$ that includes only non-zero rows is given as

$\bm{W}_{i\alpha,j\beta,k\gamma}=\frac{1}{\sqrt{2}}\begin{bmatrix}1&1&1&1&1\\\
-1&0&0&0&0\\\ 0&-1&0&0&0\\\ 0&0&-1&0&0\\\ 0&0&0&-1&0\\\ 0&0&0&0&-1\\\
\end{bmatrix}.$ (36)

Therefore, the column vector $\bm{w}_{i\alpha,j\beta,k\gamma}^{\perp}$, which
is orthonormal to all column vectors of $\bm{W}_{i\alpha,j\beta,k\gamma}$, can
be considered as an eigenvector of $\bm{P}_{w}$. It is obtained as

$\bm{w}_{i\alpha,j\beta,k\gamma}^{\perp}=\frac{1}{\sqrt{6}}[1\>1\>1\>1\>1\>1]^{\top},$
(37)

which is expressed in the submatrix representation by omitting zero elements.
If $i\alpha=j\beta$ and $i\alpha\neq k\gamma$, the permutation symmetry rule
for $(i\alpha,i\alpha,k\gamma)$ is given by Eq. (II.3). The two equations in
Eq. (II.3) are also represented as
$\bm{W}_{i\alpha,i\alpha,k\gamma}^{\top}\bm{\theta}=0$, where the submatrix of
$\bm{W}_{i\alpha,i\alpha,k\gamma}$ that includes only non-zero rows is
expressed as

$\bm{W}_{i\alpha,i\alpha,k\gamma}=\frac{1}{\sqrt{2}}\begin{bmatrix}1&1\\\
-1&0\\\ 0&-1\\\ \end{bmatrix}.$ (38)

Therefore, the vector orthonormal to all column vectors of
$\bm{W}_{i\alpha,j\beta,k\gamma}$, $\bm{w}_{i\alpha,i\alpha,k\gamma}^{\perp}$,
is given as

$\bm{w}_{i\alpha,i\alpha,k\gamma}^{\perp}=\frac{1}{\sqrt{3}}[1\>1\>1]^{\top}.$
(39)

If $i\alpha=j\beta$ and $i\alpha=k\gamma$, there exist no permutation symmetry
rules. Therefore, vector $\bm{w}_{i\alpha,i\alpha,i\alpha}^{\perp}$ is defined
as the vector that has only a non-zero element of one for column
$(i,\alpha,i,\alpha,i,\alpha)$. Thus, the eigenvectors of $\bm{P}_{w}$,
$\bm{B}_{w}$, is written as

$\bm{B}_{w}=[\cdots,\bm{w}_{i\alpha,j\beta,k\gamma}^{\perp},\cdots]^{\top},$
(40)

where vectors $\bm{w}_{i\alpha,j\beta,k\gamma}^{\perp}$ are also mutually
orthonormal. The size of $\bm{B}_{w}$ is $(27N^{3},n_{w})$, where $n_{w}\simeq
27N^{3}/6$. The number of non-zero elements is $27N^{3}$ in $\bm{B}_{w}$,
which is much smaller than the total number of elements in $\bm{B}_{w}$.

###  III.5 Effective use of lattice translation group

Here, we are introducing the projection matrix for the identity irreducible
representation of the lattice translation group, which is a normal subgroup of
the space group. Although the projection matrix for the lattice translation
group is not mandatory for describing the full projection matrix, it can help
in compressing projection matrices that are made up of the full projection
matrix. The projection matrix for the lattice translation group can be
formulated in a similar way to the derivation of the projection matrix for the
permutation symmetry rules.

The force constants remain constant for operations in lattice translation
group $\mathcal{T}=\\{t_{1},t_{2},\cdots\\}$, which are also included in group
$\mathcal{G}$. We identify representative atoms using indices
$\\{i^{\prime}\\}$ and express the translation of atom $i$ as $t_{n}(i)$. For
the second-order force constants of paired atoms, which include a
representative atom $i^{\prime}$, the invariant property is expressed as

$\Theta_{i^{\prime}\alpha,j\beta}-\Theta_{t_{n}(i^{\prime})\alpha,t_{n}(j)\beta}=0,$
(41)

where operation $t_{n}$ is not the identity. Therefore, $|\mathcal{T}|-1$
equations for each $(i^{\prime},\alpha,j,\beta)$ are represented in a matrix
form of $\bm{T}_{i^{\prime}\alpha,j\beta}$ as
$\bm{T}_{i^{\prime}\alpha,j\beta}^{\top}\bm{\theta}=0$. The submatrix of
$\bm{T}_{i^{\prime}\alpha,j\beta}$ composed only of non-zero rows is expressed
as

$\bm{T}_{i^{\prime},\alpha,j\beta}=\frac{1}{\sqrt{2}}\begin{bmatrix}1&1&1&1&\\\
-1&0&0&0&\\\ 0&-1&0&0&\cdots\\\ 0&0&-1&0&\\\ 0&0&0&-1&\\\ &\vdots&&&\ddots\\\
\end{bmatrix},$ (42)

where the submatrix is of size $(|\mathcal{T}|,|\mathcal{T}|-1)$. The vector
orthonormal to all column vectors of $\bm{T}_{i^{\prime}\alpha,j\beta}$ can be
an eigenvector of $\bm{P}_{\rm LT}$ given in the submatrix representation as

$\bm{t}_{i^{\prime}\alpha,j\beta}^{\perp}=\frac{1}{\sqrt{|\mathcal{T}|}}[1\>1\>1\>\cdots]^{\top}.$
(43)

The eigenvectors of $\bm{P}_{\rm LT}$, $\bm{B}_{\rm LT}$, is then described as

$\bm{B}_{\rm
LT}=[\cdots,\bm{t}_{i^{\prime}\alpha,j\beta}^{\perp},\cdots]^{\top},$ (44)

where vectors $\\{\bm{t}_{i^{\prime}\alpha,j\beta}^{\perp}\\}$ are also
mutually orthonormal. The size of $\bm{B}_{\rm LT}$ is
$(9N^{2},9N^{2}/|\mathcal{T}|)$. The number of non-zero elements is $9N^{2}$
in $\bm{B}_{\rm LT}$, hence, $\bm{B}_{\rm LT}$ is a sparse matrix.

The third-order force constants for triplet atoms, including a representative
atom $i^{\prime}$, are also invariant with respect to operations in
$\mathcal{T}$. This is described as

$\Theta_{i^{\prime}\alpha,j\beta,k\gamma}-\Theta_{t_{n}(i^{\prime})\alpha,t_{n}(j)\beta,t_{n}(k)\gamma}=0,$
(45)

where operation $t_{n}$ is not the identity. Similar to the second-order force
constants, $|\mathcal{T}|-1$ equations for each
$(i^{\prime},\alpha,j,\beta,k,\gamma)$ are represented in a matrix form of
$\bm{T}_{i^{\prime}\alpha,j\beta,k\gamma}$ as
$\bm{T}_{i^{\prime}\alpha,j\beta,k\gamma}^{\top}\bm{\theta}^{(\rm FC3)}=0$.
Therefore, the submatrix forms of $\bm{T}_{i^{\prime}\alpha,j\beta,k\gamma}$
and $\bm{t}_{i^{\prime}\alpha,j\beta,k\gamma}^{\perp}$ are obtained by
replacing indices ($i^{\prime}\alpha,j\beta$) with indices
($i^{\prime}\alpha,j\beta,k\gamma$) in Eqs. (42) and (43). Thus, the
eigenvectors of $\bm{P}_{\rm LT}$ is written as

$\bm{B}_{\rm
LT}=[\cdots,\bm{t}_{i^{\prime}\alpha,j\beta,k\gamma}^{\perp},\cdots]^{\top},$
(46)

where vectors $\\{\bm{t}_{i^{\prime}\alpha,j\beta,k\gamma}^{\perp}\\}$ are
also mutually orthonormal. The size of $\bm{B}_{\rm LT}$ is
$(27N^{3},27N^{3}/|\mathcal{T}|)$. The number of non-zero elements is
$27N^{3}$ in $\bm{B}_{\rm LT}$, which is much smaller than the total number of
elements in $\bm{B}_{\rm LT}$.

The projection matrix is obtained from the eigenvectors as $\bm{P}_{\rm
LT}=\bm{B}_{\rm LT}\bm{B}_{\rm LT}^{\top}$. The projection matrix has the
relationship with $\bm{P}_{\rm spg}$ of

$\bm{P}_{\rm LT}\bm{P}_{\rm spg}=\bm{P}_{\rm spg}\bm{P}_{\rm LT}=\bm{P}_{\rm
spg},$ (47)

because $t_{n}\mathcal{G}=\mathcal{G}$ for any operation $t_{n}\in\mathcal{T}$
that is also included in group $\mathcal{G}$, $t_{n}\in\mathcal{G}$.
Therefore, the projection matrix for the lattice translation group leaves
$\bm{P}_{\rm spg}$ invariant.

###  III.6 Properties of projection matrices

The commutation relations between projection matrices are significant in
building force-constant basis sets in an efficient manner. These relations
include

$[\bm{P}_{\rm spg},\bm{P}_{v}]=0,\>[\bm{P}_{\rm
spg},\bm{P}_{w}]=0,\>[\bm{P}_{v},\bm{P}_{w}]\neq 0.$ (48)

The commutation relations on the projection matrix for the lattice translation
group are also useful in compressing projection matrices in addition to the
relationship of Eq. (47). They are given as

$[\bm{P}_{\rm LT},\bm{P}_{v}]=0,\>[\bm{P}_{\rm LT},\bm{P}_{w}]=0.$ (49)

In addition, these projection matrices are idempotent and orthogonal, which
are represented by

$\bm{P}_{\rm spg}^{2}=\bm{P}_{\rm
spg},\>\bm{P}_{v}^{2}=\bm{P}_{v},\>\bm{P}_{w}^{2}=\bm{P}_{w},\>\bm{P}_{\rm
LT}^{2}=\bm{P}_{\rm LT}$ (50)

and

$\bm{P}_{\rm spg}^{\top}=\bm{P}_{\rm
spg},\>\bm{P}_{v}^{\top}=\bm{P}_{v},\>\bm{P}_{w}^{\top}=\bm{P}_{w},\>\bm{P}_{\rm
LT}^{\top}=\bm{P}_{\rm LT},$ (51)

respectively.

###  III.7 Solving eigenvalue problems for projection matrices

The techniques mentioned earlier can be used to compress projection matrices.
However, even in their compressed form, these matrices can still be quite
extensive. Nevertheless, since they are often sparse matrices, their block-
diagonal structures can be useful in solving the eigenvalue problems
associated with them in an efficient manner. In this context, we are
considering an orthogonal projection matrix that is represented in a block-
diagonal form, which is transformed by row and column permutations. The block-
diagonal form of the orthogonal projection matrix is written as

$\bm{M}^{\top}\bm{P}\bm{M}=\begin{bmatrix}\bm{I}&\bm{0}&\bm{0}\\\
\bm{0}&\bm{J}&\bm{0}\\\ \bm{0}&\bm{0}&\bm{K}\\\ \end{bmatrix},$ (52)

where $\bm{M}$ denotes the permutation matrix for column. Submatrices
$\bm{I}$, $\bm{J}$, and $\bm{K}$ must all be orthogonal projection matrices so
that $\bm{P}$ is the orthogonal projection matrix. When we solve eigenvalue
problems for the submatrices and obtain eigenvectors of $\bm{I}$, $\bm{J}$,
and $\bm{K}$ as

$\bm{I}\bm{b_{I}}=\bm{b_{I}},\>\bm{J}\bm{b_{J}}=\bm{b_{J}},\>\bm{K}\bm{b_{K}}=\bm{b_{K}},$
(53)

the vectors given by

$\bm{b}=\begin{bmatrix}\bm{b_{I}}\\\ \bm{0}\\\ \bm{0}\\\
\end{bmatrix},\begin{bmatrix}\bm{0}\\\ \bm{b_{J}}\\\ \bm{0}\\\
\end{bmatrix},\begin{bmatrix}\bm{0}\\\ \bm{0}\\\ \bm{b_{K}}\\\ \end{bmatrix}$
(54)

can be adopted as eigenvectors of $\bm{M}^{\top}\bm{P}\bm{M}$. They are
mutually orthonormal, and the trace of $\bm{M}^{\top}\bm{P}\bm{M}$ is given as

$\operatorname{tr}(\bm{M}^{\top}\bm{P}\bm{M})=\operatorname{tr}(\bm{I})+\operatorname{tr}(\bm{J})+\operatorname{tr}(\bm{K}),$
(55)

which corresponds to the number of eigenvectors of
$\bm{M}^{\top}\bm{P}\bm{M}$. Therefore, the set of eigenvectors $\\{\bm{b}\\}$
built from $\\{\bm{b_{I}}\\}$, $\\{\bm{b_{J}}\\}$, and $\\{\bm{b_{K}}\\}$ is
complete. The eigenvectors of $\bm{P}$ are finally calculated as
$\bm{M}\bm{b}$.

We follow a graph-based approach to determine the permutation matrix $\bm{M}$.
To achieve this, we treat the orthogonal projection matrix $\bm{P}$ as an
undirected graph where each row and column of the matrix corresponds to a
vertex. The vertices $s$ and $t$ are considered adjacent only when the element
of $\bm{P}$ for row $s$ and column $t$ is non-zero. Next, we employ an
efficient graph algorithm to identify the connected components Sedgewick
(2001) in the graph representation of the projection matrix. The algorithm
assigns columns and rows of the projection matrix to one of the connected
components, which represent block structures in the projection matrix.
Finally, we construct the permutation matrix $\bm{M}$ such that the columns in
the same connected component are arranged adjacent to each other in the
reordered projection matrix $\bm{M}^{\top}\bm{P}\bm{M}$.

###  III.8 Procedure to calculate eigenvectors of full projection matrix

The eigenvectors of the approximated projection matrix $\tilde{\bm{P}}$ given
by Eq. (27) is calculated by following the equations shown below. We start by
rewriting matrix $\tilde{\bm{P}}$ using the commutation relations and the
idempotent properties of projection matrices as

$\displaystyle\tilde{\bm{P}}$ $\displaystyle=$ $\displaystyle\bm{P}_{\rm
spg}\bm{P}_{w}\bm{P}_{v}\bm{P}_{w}$ (56) $\displaystyle=$
$\displaystyle\bm{P}_{\rm spg}\bm{P}_{w}\bm{P}_{v}\bm{P}_{\rm spg}\bm{P}_{w}.$

The eigenvalue decomposition of the matrix product $\bm{P}_{\rm
spg}\bm{P}_{w}$ is calculated as

$\displaystyle\bm{P}_{\rm spg}\bm{P}_{w}$ $\displaystyle=$
$\displaystyle\bm{P}_{\rm LT}\bm{P}_{\rm spg}\bm{P}_{\rm
LT}\bm{P}_{w}\bm{P}_{\rm LT}$ (57) $\displaystyle=$ $\displaystyle\bm{B}_{\rm
LT}\bm{B}_{\rm LT}^{\top}\bm{P}_{\rm spg}\bm{B}_{\rm LT}\bm{B}_{\rm
LT}^{\top}\bm{B}_{w}\bm{B}_{w}^{\top}\bm{B}_{\rm LT}\bm{B}_{\rm
LT}^{\top}\>\>\>(\bm{P}_{\rm LT}=\bm{B}_{\rm LT}\bm{B}_{\rm
LT}^{\top},\bm{P}_{w}=\bm{B}_{w}\bm{B}_{w}^{\top})$ $\displaystyle=$
$\displaystyle\bm{B}_{\rm LT}\bm{P}_{\rm spg\downarrow
LT}\bm{P}_{w\downarrow{\rm LT}}\bm{B}_{\rm LT}^{\top}\>\>\>\text{(Eqs.
(\ref{symfc:matrix-comp-definition}) and (\ref{symfc:matrix-comp-
definition2}))}$ $\displaystyle=$ $\displaystyle\bm{B}_{\rm LT}\bm{P}_{\rm
spg\downarrow LT}\bm{B}_{w\downarrow{\rm LT}}\bm{B}_{w\downarrow{\rm
LT}}^{\top}\bm{B}_{\rm LT}^{\top}\>\>\>(\bm{P}_{w\downarrow{\rm
LT}}=\bm{B}_{w\downarrow{\rm LT}}\bm{B}_{w\downarrow{\rm LT}}^{\top})$
$\displaystyle=$ $\displaystyle\bm{B}_{\rm LT}(\bm{B}_{w\downarrow{\rm
LT}}\bm{B}_{w\downarrow{\rm LT}}^{\top})\bm{P}_{\rm spg\downarrow
LT}\bm{B}_{w\downarrow{\rm LT}}\bm{B}_{w\downarrow{\rm LT}}^{\top}\bm{B}_{\rm
LT}^{\top}\>\>\>$ $\displaystyle=$ $\displaystyle\bm{B}_{\rm
LT}\bm{B}_{w\downarrow{\rm LT}}\bm{P}_{({\rm spg\downarrow
LT})\downarrow(w\downarrow{\rm LT})}\bm{B}_{w\downarrow{\rm
LT}}^{\top}\bm{B}_{\rm LT}^{\top}\>\>\>\text{(Eqs. (\ref{symfc:matrix-comp-
definition}) and (\ref{symfc:matrix-comp-definition2}))}$ $\displaystyle=$
$\displaystyle\bm{B}_{\rm LT}\bm{B}_{w\downarrow{\rm LT}}\bm{B}_{({\rm
spg\downarrow LT})\downarrow(w\downarrow{\rm LT})}\bm{B}_{({\rm spg\downarrow
LT})\downarrow(w\downarrow{\rm LT})}^{\top}\bm{B}_{w\downarrow{\rm
LT}}^{\top}\bm{B}_{\rm LT}^{\top}$
$\displaystyle\>\>\>\>\>\>\>\>\>\>(\bm{P}_{({\rm spg\downarrow
LT})\downarrow(w\downarrow{\rm LT})}=\bm{B}_{({\rm spg\downarrow
LT})\downarrow(w\downarrow{\rm LT})}\bm{B}_{({\rm spg\downarrow
LT})\downarrow(w\downarrow{\rm LT})}^{\top})$ $\displaystyle=$
$\displaystyle\bm{B}_{{\rm spg}\cap w}\bm{B}_{{\rm spg}\cap w}^{\top},$

where

$\bm{B}_{{\rm spg}\cap w}=\bm{B}_{\rm LT}\bm{B}_{w\downarrow{\rm
LT}}\bm{B}_{({\rm spg\downarrow LT})\downarrow(w\downarrow{\rm LT})}.$ (58)

Then, the eigenvalue decomposition of matrix $\tilde{\bm{P}}$ is computed by

$\displaystyle\tilde{\bm{P}}$ $\displaystyle=$ $\displaystyle\bm{P}_{\rm
spg}\bm{P}_{w}\bm{P}_{v}\bm{P}_{\rm spg}\bm{P}_{w}$ (59) $\displaystyle=$
$\displaystyle\bm{B}_{{\rm spg}\cap w}\bm{B}_{{\rm spg}\cap
w}^{\top}\bm{P}_{v}\bm{B}_{{\rm spg}\cap w}\bm{B}_{{\rm spg}\cap w}^{\top}$
$\displaystyle=$ $\displaystyle\bm{B}_{{\rm spg}\cap
w}\bm{P}_{v\downarrow({\rm spg}\cap w)}\bm{B}_{{\rm spg}\cap w}^{\top}$
$\displaystyle=$ $\displaystyle\bm{B}_{{\rm spg}\cap
w}\bm{B}_{v\downarrow({\rm spg}\cap w)}\bm{B}_{v\downarrow({\rm spg}\cap
w)}^{\top}\bm{B}_{{\rm spg}\cap w}^{\top}$ $\displaystyle=$
$\displaystyle\tilde{\bm{B}}\tilde{\bm{B}}^{\top},$

where the basis set for the force constants are

$\tilde{\bm{B}}=\bm{B}_{{\rm spg}\cap w}\bm{B}_{v\downarrow({\rm spg}\cap
w)}.$ (60)

Note that eigenvectors of $\tilde{\bm{P}}$ corresponding to its eigenvalues
that are less than one are eliminated, as described above.

##  IV Force constant estimation

###  IV.1 Force-constant basis set expansion

The forces acting on atoms $\\{f_{i\alpha}^{(s)}\\}$ in a supercell structure
$s$ corresponding to atomic displacements $\\{u_{i\alpha}^{(s)}\\}$ can be
simply written as

$f_{i\alpha}^{(s)}=-\sum_{j\beta}\Theta_{i\alpha,j\beta}u_{j\beta}^{(s)}-\frac{1}{2}\sum_{j\beta}\sum_{k\gamma}\Theta_{i\alpha,j\beta,k\gamma}u_{j\beta}^{(s)}u_{k\gamma}^{(s)}+\cdots.$
(61)

The force constants are expanded in basis sets given by Eq. (60) as

$\displaystyle\Theta_{i\alpha,j\beta}$ $\displaystyle=$
$\displaystyle\sum_{b}w_{b}^{\rm(FC2)}\tilde{B}_{i\alpha,j\beta,b}^{\rm(FC2)}$
(62) $\displaystyle=$
$\displaystyle\sum_{b}w_{b}^{\rm(FC2)}\left[\sum_{c}C_{i\alpha,j\beta,c}^{\rm(FC2)}E_{c,b}^{\rm(FC2)}\right]$

and

$\displaystyle\Theta_{i\alpha,j\beta,k\gamma}$ $\displaystyle=$
$\displaystyle\sum_{b}w_{b}^{\rm(FC3)}\tilde{B}_{i\alpha,j\beta,k\gamma,b}^{\rm(FC3)}$
(63) $\displaystyle=$
$\displaystyle\sum_{b}w_{b}^{\rm(FC3)}\left[\sum_{c}C_{i\alpha,j\beta,k\gamma,c}^{\rm(FC3)}E_{c,b}^{\rm(FC3)}\right],$

where $w_{b}$ denotes the expansion coefficient for $b$-th basis. Here, the
symbols $\bm{B}_{{\rm spg}\cap w}$ and $\bm{B}_{v\downarrow({\rm spg}\cap w)}$
shown in Eq. (60) are abbreviated as $\bm{C}$ and $\bm{E}$, respectively.
Then, the forces in structure $s$ are expressed as

$\displaystyle f_{i\alpha}^{(s)}=$ $\displaystyle-$
$\displaystyle\sum_{b}w_{b}^{\rm(FC2)}\sum_{j\beta}\sum_{c}u_{j\beta}^{(s)}C_{i\alpha,j\beta,c}^{\rm(FC2)}E_{c,b}^{\rm(FC2)}$
(64) $\displaystyle-$
$\displaystyle\frac{1}{2}\sum_{b}w_{b}^{\rm(FC3)}\sum_{j\beta}\sum_{k\gamma}\sum_{c}u_{j\beta}^{(s)}u_{k\gamma}^{(s)}C_{i\alpha,j\beta,k\gamma,c}^{\rm(FC3)}E_{c,b}^{\rm(FC3)}$
$\displaystyle+$ $\displaystyle\cdots.$

Therefore, the forces are rewritten as

$\displaystyle
f_{i\alpha}^{(s)}=\sum_{b}w_{b}^{\rm(FC2)}x_{i\alpha,b}^{(s,{\rm
FC2})}+\sum_{b}w_{b}^{\rm(FC3)}x_{i\alpha,b}^{(s,{\rm FC3})}+\cdots,$ (65)

where

$\displaystyle x_{i\alpha,b}^{(s,{\rm
FC2})}=-\sum_{j\beta}\sum_{c}u_{j\beta}^{(s)}C_{i\alpha,j\beta,c}^{\rm(FC2)}E_{c,b}^{\rm(FC2)}$
(66)

and

$\displaystyle x_{i\alpha,b}^{(s,{\rm
FC3})}=-\frac{1}{2}\sum_{j\beta}\sum_{k\gamma}\sum_{c}u_{j\beta}^{(s)}u_{k\gamma}^{(s)}C_{i\alpha,j\beta,k\gamma,c}^{\rm(FC3)}E_{c,b}^{\rm(FC3)}.$
(67)

The forces in structure $s$ are also expressed in a matrix form by considering
the force constants up to the third order as

$\bm{f}^{(s)}=\bm{X}^{(s)}\bm{w},$ (68)

where

$\bm{f}^{(s)}=\begin{bmatrix}\vdots\\\ f_{i\alpha}^{(s)}\\\ \vdots\\\
\end{bmatrix},$ (69)

$\displaystyle\bm{X}^{(s)}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}\bm{X}^{(s,\rm FC2)}&\bm{X}^{(s,\rm FC3)}\\\
\end{bmatrix}$ (70) $\displaystyle=$
$\displaystyle\begin{bmatrix}\cdots&\vdots&\cdots&\vdots&\cdots\\\
\cdots&x_{i\alpha,b}^{(s,{\rm FC2})}&\cdots&x_{i\alpha,b}^{(s,{\rm
FC3})}&\cdots\\\ \cdots&\vdots&\cdots&\vdots&\cdots\\\ \end{bmatrix},$

and

$\bm{w}=\begin{bmatrix}\bm{w}^{(\rm FC2)}\\\ \bm{w}^{(\rm FC3)}\\\
\end{bmatrix}=\begin{bmatrix}\vdots\\\ w_{b}^{\rm(FC2)}\\\ \vdots\\\
w_{b}^{\rm(FC3)}\\\ \vdots\\\ \end{bmatrix}.$ (71)

###  IV.2 Force constant estimation

The expansion coefficients $\bm{w}$ are determined from the forces acting on
atoms in $N_{s}$ structures, which correspond to the displacement set, using
standard linear regression. The training dataset is expressed in a matrix form
as

$\bm{X}=\begin{bmatrix}\bm{X}^{(1)}\\\ \bm{X}^{(2)}\\\ \bm{X}^{(3)}\\\
\vdots\\\
\bm{X}^{(N_{s})}\end{bmatrix},\>\>\bm{y}=\begin{bmatrix}\bm{f}^{(1)}\\\
\bm{f}^{(2)}\\\ \bm{f}^{(3)}\\\ \vdots\\\ \bm{f}^{(N_{s})}\end{bmatrix},$ (72)

and the solution for the coefficients $\bm{w}$ is given by

$\bm{w}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}\bm{y}.$ (73)

The force constants can then be recovered from the coefficients $\bm{w}$ as

$\displaystyle\bm{\theta}^{(\rm FC2)}=\bm{C}^{(\rm FC2)}\bm{E}^{(\rm
FC2)}\bm{w}^{(\rm FC2)}$ $\displaystyle\bm{\theta}^{(\rm FC3)}=\bm{C}^{(\rm
FC3)}\bm{E}^{(\rm FC3)}\bm{w}^{(\rm FC3)}.$ (74)

However, computing matrix $\bm{X}$ can be a resource-intensive process,
especially when dealing with large displacement-force datasets and basis sets.
To optimize memory allocation and efficiently estimate expansion coefficients,
we employ a technique from the recursive method for linear regression. We
begin by considering the normal equation in linear regression,
$\bm{X}^{\top}\bm{X}\bm{w}=\bm{X}^{\top}\bm{y}$. When the dataset is divided
into $n_{q}$ subsets, $\bm{X}^{\top}\bm{X}$ and $\bm{X}^{\top}\bm{y}$ can be
calculated as

$\bm{X}^{\top}\bm{X}=\sum_{q=1}^{n_{q}}\bm{X}_{q}^{\top}\bm{X}_{q},\>\>\>\bm{X}^{\top}\bm{y}=\sum_{q=1}^{n_{q}}\bm{X}_{q}^{\top}\bm{y}_{q},$
(75)

where $\bm{X}_{q}$ and $\bm{y}_{q}$ are $q$-th subsets of $\bm{X}$ and
$\bm{y}$, respectively. This allows us to avoid storing the entire $\bm{X}$
matrix. The coefficients are obtained by solving the normal equation with fast
linear algebra libraries. When dividing the displacement-force dataset, we
adopt a procedure to compute $\bm{X}^{\top}\bm{X}$ and $\bm{X}^{\top}\bm{y}$
more efficiently, which is shown in the Appendix C.

Figure 1:  (a) Number of basis vectors that completely represent second-order
force constants in elemental crystal structures, i.e., the body-centered cubic
(bcc), face-centered cubic (fcc), hexagonal close packed (hcp), simple cubic
(sc), and diamond-type structures. Numbers of basis vectors that completely
represent third-order force constants in (b) elemental crystal structures and
(c) ionic crystal structures are also shown. The ionic crystal structures
consist of the rocksalt, zincblende, wurtzite, fluorite, rutile, cubic
perovskite (C-perovskite), and spinel structures. The black solid lines
indicate the numbers of $9N^{2}$ and $27N^{3}$ for the second order and third
order, respectively.

##  V Results and discussion

###  V.1 Basis set construction

Given a unit cell structure and a supercell size, it is possible to
efficiently calculate a basis set of second-order and third-order supercell
force constants using Eqs. (57) and (59). This basis set represents supercell
force constants that follow the permutational symmetry rules, sum rules, and
symmetric properties of the given unit cell. Figure 1 (a) shows the size of
the basis set for the second-order force constants of several crystal
structures found in elemental systems. The total number of second-order force
constant elements is 9$N^{2}$, but they are significantly reduced by imposing
constraints on the force constants. For the second order, we can obtain basis
sets at small computational costs for supercells with up to 5000 atoms, which
have approximately $10^{8}$ force constant elements.

The total number of third-order force-constant elements is 27$N^{3}$, which
can become quite large when evaluating force constants using a typical size of
supercells. However, there is an efficient procedure available that can help
reduce the force constants and create orthonormal basis sets. Figures 1 (b)
and (c) show the number of basis vectors for the third-order force constants
of crystal structures found in elemental and ionic systems, respectively. For
elemental structures, the current procedure can be used with supercells of up
to 900 atoms using a common workstation, thanks to their small primitive cells
and high symmetric properties. In the case of the ionic structures, the
applicable numbers of atoms in supercells are reduced by decreasing lattice
translations in these structures. For the well-known structures shown in Fig.
1 (c), the current procedure usually works well for supercells with up to 300
atoms. More information on the corresponding computational costs for
constructing the basis sets can be found in the Appendix A.

###  V.2 Force constant evaluation

We will present four different applications of the current approach for
estimating force constants. The first application involves using a machine
learning potential (MLP) to calculate the lattice thermal conductivities
(LTCs) in diamond Si. The second application is to calculate LTCs using the
the density functional theory (DFT) calculation in diamond Si. The third
application is to calculate LTCs using an MLP in wurtzite AgI. Finally, we use
the DFT calculation to calculate LTCs in wurtzite-type ZnS. These applications
demonstrate the importance of choosing reasonable parameters and provide
recommended settings for this approach.

####  V.2.1 Computational details

Supercell force constants are determined by applying standard linear
regression to displacement-force datasets. A displacement-force dataset
consists of $N_{s}$ supercell structures that are created by introducing
random atomic displacements of magnitude $d$ into the original supercell
structure with a given supercell size. The forces acting on atoms in these
structures are then computed using MLPs and DFT calculations. In applications
using MLPs, the forces acting on atoms in these structures are computed using
polynomial MLPs Seko  _et al._ (2019); Seko (2023) due to their negligible
computational costs in calculating force constants. The polynomial MLP
represents short-range interatomic interactions with a polynomial function of
polynomial invariants for arbitrary rotations. In applications using DFT
calculations, force calculations were performed using the plane-wave-basis
projector augmented wave method Blöchl (1994) within the
Perdew–Burke–Ernzerhof exchange-correlation functional Perdew  _et al._ (1996)
as implemented in the vasp code Kresse and Hafner (1993); Kresse and
Furthmüller (1996); Kresse and Joubert (1999).

The LTCs were calculated as the solution of the Peierls–Boltzmann equation
under the relaxation time approximation Peierls (1929, 2001); Allen and
Perebeinos (2018) using the phono3py code Togo  _et al._ (2015, 2023). To
simplify the LTC calculations and measure the precision of force constants,
only phonon-phonon scattering was considered to obtain the phonon relaxation
times. The phonon lifetimes were calculated from supercell force constants as
the reciprocals of the imaginary part of the phonon self-energy corresponding
to the bubble diagram, which were used as the relaxation times. Phonon
properties required for the LTC calculations such as phonon group velocities
and mode heat capacities were obtained from dynamical matrices derived from
second-order supercell force constants. Mesh grids for sampling the reciprocal
spaces are shown in each application.

####  V.2.2 Diamond Si (MLP)

Here, we use the same MLP as that used in the literature Fujii and Seko
(2022), where LTCs at grain boundaries in Si are studied using molecular
dynamics simulations. This MLP is also available in the Polynomial Machine
Learning Potential Repository Mac , and it is estimated from a dataset with
various crystal structures and their derivatives where random atomic
displacements and random cell shape changes are introduced. In addition, the
reciprocal spaces for the LTC calculations were sampled by the $19\times
19\times 19$ mesh.

Figure 2:  (a) Relationship between the LTC and the number of supercells with
$d=0.001$, $0.003$, and $0.01$ Å in diamond Si. The minimum numbers of
supercells required for constructing a full-rank predictor matrix are 5, 14,
and 33 for the $2\times 2\times 2$, $3\times 3\times 3$, and $4\times 4\times
4$ supercell expansions, respectively. Only LTC values for the number of
supercells up to 100 are shown. (b) Relationship between the LTC and the
number of supercells with $d=0.03$, $0.05$, and $0.1$ Å in diamond Si.

The current approach based on projectors allows us to efficiently estimate
force constants and LTCs for supercells with up to 512 atoms. Figure 2 shows
the relationship between the LTC, the displacement magnitude, and the number
of supercell structures in the displacement-force dataset. The complete basis-
set sizes for the $2\times 2\times 2$, $3\times 3\times 3$, and $4\times
4\times 4$ supercell expansions are 777, 8800, and 49301, respectively.
Therefore, the minimum numbers of supercells needed to build a full-rank
predictor matrix $\bm{X}$ are 5, 14, and 33, respectively. In Fig. 2 (a), it
is evident that even a small number of supercells can produce precise LTCs
when small displacement values are considered. Specifically, 5 and 15
supercell structures generated from random displacements with $d=0.001$ Å are
sufficient to derive precise LTCs for the $2\times 2\times 2$ and $3\times
3\times 3$ expansions, respectively. Conversely, the minimal structures are
not able to precisely predict the LTC when larger displacements are employed.
The LTC converges slowly as the number of supercells increases as can be seen
in Fig. 2 (b).

Figure 3:  Relationship between the LTC value and the magnitude of atomic
displacements in diamond Si with the supercell size of $2\times 2\times 2$.
The LTC values are calculated using displacement-force datasets composed of
10, 100, 1000, and 10000 supercell structures.

Figure 3 shows the relationship between the LTC value and the magnitude of
atomic displacements. When using a quite small value of displacements, such as
$d=10^{-6}$ Å, the forces exhibit quite small values, which leads to numerical
noise becoming significant enough to contribute to force constants. As a
result, many supercells are required to obtain a precise LTC. For accurate
evaluation of third-order force constants without any higher-order
contributions, it is better to choose the magnitude of atomic displacements
within the range of $10^{-5}$ to $10^{-3}$ Å.

Similarly, when using atomic displacements larger than 0.01 Å, a large number
of supercells are also necessary to obtain the converged LTC. This is due to
the truncation errors of higher-order contributions, which make it
computationally challenging to calculate the converged LTC. It is worth noting
that we have utilized a constant magnitude of atomic displacements to generate
each displacement-force dataset in this study. However, the current discussion
is relevant to displacement-force datasets that consist of large and varying
magnitudes of atomic displacements. Such datasets are found in estimating
temperature-dependent force constants relevant to self-consistent phonons,
which effectively contain higher-order contributions of force constants.
Therefore, if self-consistent phonons are evaluated at a high temperature, we
will need to use a large number of supercell structures for the displacement-
force dataset.

Figure 4:  Dependence of the condition number of $\bm{X}^{\top}\bm{X}$ for the
third-order force constants on the number of supercells in diamond Si. The
supercell size of $2\times 2\times 2$ and the magnitude of atomic
displacements of $d=0.001$ Å are used. The condition number is almost
independent of the magnitude of atomic displacements, and similar values of
the condition number can be found when using $d=0.1$ Å.

The imprecise prediction of LTC is caused by errors in the forces, which are
strongly related to linearly dependent relationships between predictor
variables in matrix $\bm{X}$. This is known as multicollinearity, and it can
make the regression coefficients sensitive to force errors. In standard linear
regression, the condition number of matrix $\bm{X}^{\top}\bm{X}$ is typically
used as a measure of multicollinearity. The condition number is calculated by
dividing the maximum eigenvalue of $\bm{X}^{\top}\bm{X}$ by the minimum
eigenvalue. If the condition number is significantly large, multicollinearity
is likely to occur in regression.

Figure 4 shows the dependence of the condition number of matrix
$\bm{X}^{\top}\bm{X}$ for the third-order force constants with the $2\times
2\times 2$ supercell on the number of supercell structures. The condition
number is high when using only a small number of supercells, which is close to
the minimum number of supercells required for constructing a full-rank
predictor matrix. This means that the prediction of LTC values is sensitive to
force errors, resulting in imprecise predictions of LTC values. On the other
hand, the condition number is lower at larger numbers of supercells and
gradually decreases as the number of supercells increases. This behavior is
consistent with the convergence patterns in LTC. Thus, by using both an
orthonormal basis set for force constants and random displacements for force
calculations, it is possible to create optimal conditions for estimating force
constants through standard linear regression.

####  V.2.3 Diamond Si (DFT)

Figure 5:  (a) Relationship between the LTC and the number of supercells with
$d=0.003$, $0.01$, $0.03$, and $0.1$ Å in diamond Si. Displacement-force
datasets are developed using only the DFT calculation. The minimum number of
supercells required for constructing a full-rank predictor matrix is five for
the $2\times 2\times 2$ supercell expansion. Only LTC values for the number of
supercells up to 50 are shown. (b) Relationship between the LTC and the number
of supercells with $d=0.003$, $0.01$, $0.03$, and $0.1$ Å in diamond Si. LTC
values for the number of supercells up to 1000 are shown.

We will demonstrate an application to diamond Si, where displacement-force
datasets are developed using only DFT calculations. In this application, we
computed the forces acting on atoms for supercell structures using DFT
calculations, and used them to estimate force constants. For the DFT
calculations, we constructed supercells using the $2\times 2\times 2$
expansion of the unit cell. The cutoff energy was set to 300 eV, and the total
energies converged to less than 10-3 meV/supercell. We considered the
magnitudes of atomic displacements as $d=0.003$, $0.01$, $0.03$, and $0.1$ Å.

In Fig. 5, we present the relationship between the LTC and the number of
supercells in the displacement-force dataset. The LTC values obtained from the
DFT calculation are slightly larger than those obtained from the polynomial
MLP, as demonstrated in Ref. Fujii and Seko, 2022. When atomic displacements
are larger than 0.01 Å and the number of supercells is small, the LTC values
are underestimated due to errors in the forces resulting from the truncation
of higher-order contributions. Therefore, the minimal supercell structures
fail to predict the LTC precisely.

####  V.2.4 Wurtzite AgI (MLP)

To calculate LTCs in wurtzite AgI, we use a polynomial MLP which was developed
from a displacement-force DFT dataset using an on-the-fly approach Togo and
Seko (2024). This dataset was made up of 300 supercell structures with the
magnitude of atomic displacements of 0.03 Å. We assess the forces acting on
atoms for $N_{s}$ structures using the polynomial MLP, where residual forces
for the calculated equilibrium atomic positions are eliminated to stabilize
the LTC calculations with small atomic displacements. Additionally, we sample
the reciprocal spaces for the LTC calculations using the $19\times 19\times
10$ mesh.

Figure 6:  Convergence behavior of LTC with respect to the number of
supercells in wurtzite AgI. The supercell size is $3\times 3\times 2$, and the
number of basis vectors for the third order force constants is 7752. The
average values of diagonal elements in LTC tensors are shown.

Figure 6 shows the relationship between the LTC and the number of supercells
in wurtzite AgI. The supercells are built using the $3\times 3\times 2$
expansion of the unit cell. To construct a full-rank predictor matrix
$\bm{X}$, a minimum of 36 supercell structures is needed, given that the
number of basis vectors for the third-order force constants is 7752. The
convergence behavior of the LTC is determined by the magnitude of atomic
displacements. For instance, a small number of supercell structures is
sufficient to obtain a converged value of LTC for small atomic displacements,
similar to the case of diamond Si. However, for large displacements, a larger
number of structures are required to obtain a precise LTC value, where higher-
order contributions are effectively included.

####  V.2.5 Wurtzite ZnS (DFT)

Figure 7:  Convergence behavior of LTC with respect to the number of
supercells in the wurtzite-type ZnS. Displacement-force datasets are developed
using only the DFT calculation. The supercell size is $3\times 3\times 2$, and
the minimum number of supercells required for constructing a full-rank
predictor matrix is 36. The average values of diagonal elements in LTC tensors
are shown. The LTC values calculated using the forces from which the residual
forces are not subtracted are represented by the closed circles, while the LTC
values calculated using the forces where the residual forces are eliminated
are represented by the open circles.

We are demonstrating LTC calculations using the DFT calculation in wurzite-
type ZnS. Displacement-force datasets were created using only DFT
calculations. To estimate force constants, we utilized DFT calculations to
compute the forces acting on atoms in supercell structures. For these
calculations, we considered supercells constructed using the $3\times 3\times
2$ expansion of the unit cell. The cutoff energy was set to 500 eV, and the
total energies converged to less than 0.001 meV/supercell. The atomic
positions in the unit cell were optimized until the residual forces were less
than $0.0015$ eV/Å. To examine the atomic displacement dependence of the LTC,
we used magnitudes of atomic displacements of $d=0.003$, $0.01$, $0.03$, and
$0.1$ Å.

Figure 7 shows the relationship between the LTC, the number of supercells, and
the magnitude of atomic displacements in the wurtzite-type ZnS. The DFT
calculation is used to determine the equilibrium structure of the wurtzite
structure, but residual forces still remain in the equilibrium structure
within a given convergence criterion. These residual forces become a source of
errors in the forces and can have a significant impact on the prediction of
the LTC, particularly when the magnitudes of atomic displacements and forces
are small. In the case of using $d=0.003$ and $0.01$ Å, the LTC values are
sensitive to whether or not the residual forces are subtracted from the DFT
values of the forces in supercell structures, as seen in Fig. 7. On the other
hand, the magnitudes of the forces in supercell structures with $d=0.03$ and
$0.1$ Å are large enough that the LTC values are independent of the
subtraction of residual forces.

From the results in Fig. 7, it is clear that subtracting residual forces is
essential for estimating force constants and LTCs precisely and robustly.
However, when using $d=0.003$ Å, the converged LTC value is slightly different
from those for $d=0.01$ and $0.03$ Å. This difference may arise from the fact
that the residual forces and other numerical errors in the DFT calculation are
too large to estimate the LTC precisely. On the other hand, when using $d=0.1$
Å, higher-order contributions are significant, and their truncation causes the
slow convergence of the LTC with respect to the number of supercells.
Therefore, it is recommended that the magnitude of atomic displacements is set
between $d=0.01$ and $0.03$ Å in this case.

In comparison to the diamond structure, a greater number of supercell
structures are needed due to the increase of basis vectors for the force
constants. In such a case, it is beneficial to adopt a more efficient approach
for determining the force constants Togo and Seko (2024). In this approach,
the current projector-based procedure is combined with the use of on-the-fly
MLPs, which are developed from the displacement-force dataset. This approach
can significantly reduce the number of supercell structures required for the
DFT calculation. At the same time, the force constants and LTCs can be
predicted as accurately as those calculated from DFT calculations for a much
larger number of structures.

##  VI Conclusion

In this study, we have proposed a projector-based approach to efficiently
build a complete orthonormal basis set for the second-order and third-order
force constants. We have introduced several compression techniques that enable
the reduction of the basis set in an efficient manner. These techniques
compress a projection matrix using eigenvectors of another projection matrix
onto subspaces defined by partial requirements for the force constants. Our
efficient procedures for the basis-set construction and force constant
estimation from a displacement-force dataset are implemented in the symfc code
Seko and Togo . Additionally, extending this approach to higher-order force
constants may be straightforward.

The obtained basis set satisfies all the necessary requirements for the force
constants, including the permutation symmetry rules, sum rules, and symmetric
properties required by the space group operations of the target structure.
Therefore, the complete orthonormal basis set helps efficiently obtain precise
force constants that exactly satisfy constraints on the force constants. As a
result, properties related to the force constants, such as LTCs, can be
precisely and efficiently calculated using the basis set for the force
constants, as demonstrated in the applications of this study.

###### Acknowledgements.

This work was supported by JSPS KAKENHI Grant Numbers JP22H01756, and
JP19H05787, and JP21K04632.

##  Appendix A Computational efficiency

Figure 8:  Elapsed time required for computing basis sets for elemental and
ionic structures. The elapsed times are measured using Intel(R) Xeon(R) CPU
E5-2683 v4.

In this section, we present benchmark results that demonstrate the
computational efficiency of generating force-constant basis sets and computing
force constants from displacement-force datasets. Figure 8 shows the elapsed
time required for computing the basis set of second-order and third-order
supercell force constants for elemental and ionic structures. The
computational cost mainly depends on the number of atoms in the supercell
because the numbers of second-order and third-order force-constant elements
increase as $9N^{2}$ and $27N^{3}$, respectively. Moreover, the computational
cost also varies depending on the symmetry of the target structure. Structures
with lower symmetry tend to require a higher computational cost than those
with higher symmetry, since their symmetric properties can efficiently reduce
basis sets. This is ascribed to the fact that it is more computationally
demanding to solve the compressed eigenvalue problem given by Eq. (59) in the
former.

Figure 9:  (a) Elapsed time and peak memory usage in estimating force
constants with different subset sizes in diamond Si. The number of supercells,
the supercell size, and the magnitude of atomic displacements are 3000,
$3\times 3\times 3$ (216 atoms), and $d=0.001$ Å, respectively. (b) Elapsed
time and peak memory usage in estimating force constants with a fixed subset
size of 100 and different numbers of supercells in diamond Si. The subset
size, the supercell size, and the magnitude of atomic displacements are 3000,
$3\times 3\times 3$ (216 atoms), and $d=0.001$ Å, respectively. (c) Elapsed
time and peak memory usage in estimating force constants with a fixed subset
size of 100 and a fixed number of supercells of 3000. The magnitude of atomic
displacements is set to $d=0.001$ Å. The elapsed times are measured using
Intel(R) Xeon(R) CPU E5-2683 v4.

Then, we have conducted a benchmark to estimate force constants from
displacement-force datasets in diamond Si. Figure 9 (a) shows the elapsed time
and peak memory usage in estimating force constants with different subset
sizes. For this benchmark, we use a dataset comprising 3000 structures with a
supercell size of $3\times 3\times 3$ in diamond Si. As shown in Fig. 9 (a),
the computational cost does not depend on the subset size of up to 400
structures, while the required memory size is almost proportional to the
subset size. Therefore, we suggest that the subset size should not exceed 100.

Figure 9 (b) shows the elapsed time and peak memory usage in estimating force
constants with a fixed subset size of 100, but with different numbers of
supercells. The computational cost is proportional to the number of supercells
in the dataset, while the required memory size remains almost constant,
irrespective of the number of supercells. This result indicates that the
subset division of $\bm{X}^{\top}\bm{X}$ is effective in reducing memory
consumption.

Finally, in Fig. 9 (c), we present the elapsed time and peak memory usage in
estimating force constants with a fixed subset size of 100 and a fixed number
of supercells of 3000. We use three supercell sizes of $2\times 2\times 2$,
$3\times 3\times 3$, and $4\times 4\times 4$ for the benchmark. The results
indicate that the computational cost and the required memory size strongly
depend on the number of atoms in the supercell.

##  Appendix B Commutation relations of projection matrices

In this section, we present a proof of the commutation relations between
projection matrices, $[\bm{P}_{\rm spg},\bm{P}_{v}]=0$ and $[\bm{P}_{\rm
spg},\bm{P}_{w}]=0$. Although we only provide a proof for the second-order
force constants, it can easily be extended to the third-order force constants.
Suppose we have two orthogonal projection matrices, $\bm{P}_{1}$ and
$\bm{P}_{2}$. Therefore, the transpose of the product of $\bm{P}_{1}$ and
$\bm{P}_{2}$ is given by
$(\bm{P}_{1}\bm{P}_{2})^{\top}=\bm{P}_{2}^{\top}\bm{P}_{1}^{\top}=\bm{P}_{2}\bm{P}_{1}$.
If $\bm{P}_{1}\bm{P}_{2}$ is symmetric, i.e.,
$\bm{P}_{1}\bm{P}_{2}=(\bm{P}_{1}\bm{P}_{2})^{\top}$, then $\bm{P}_{1}$ and
$\bm{P}_{2}$ commute. In the following, we will demonstrate that the product
of two projection matrices is symmetric.

(1) Proof of $[\bm{P}_{\rm spg},\bm{P}_{v}]=0$. Firstly, we will show that
$\bm{P}_{v}\bm{P}_{\rm spg}$ is symmetric, $\bm{P}_{v}\bm{P}_{\rm
spg}=(\bm{P}_{v}\bm{P}_{\rm spg})^{\top}$. Since the projection matrix
$\bm{P}_{v}$ is given by $\bm{P}_{v}=\bm{I}-\bm{C}_{v}\bm{C}_{v}^{\top}$ as
shown in Eq. (22), the elements of
$\bm{P}_{v}^{c}=\bm{C}_{v}\bm{C}_{v}^{\top}$ are expressed as

$P^{c}_{v,(i\alpha
j\beta,i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime})}=\frac{1}{N}\delta_{ii^{\prime}}\delta_{\alpha\alpha^{\prime}}\delta_{\beta\beta^{\prime}},$
(76)

where $\delta_{ii^{\prime}}$ denotes the Kronecker delta. Therefore, when
denoting the product of $\bm{P}_{v}^{c}$ and $\bm{P}_{\rm spg}$ as
$\bm{P}^{\prime}=\bm{P}_{v}^{c}\bm{P}_{\rm spg}$, its elements are written as

$\displaystyle P^{\prime}_{(i\alpha
j\beta,i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime})}$ $\displaystyle=$
$\displaystyle\sum_{k\gamma l\delta}P^{c}_{v,(i\alpha j\beta,k\gamma
l\delta)}P_{{\rm spg}{(k\gamma
l\delta,i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime})}}$ (77)
$\displaystyle=$ $\displaystyle\frac{1}{N}\sum_{l}P_{{\rm spg}{(i\alpha
l\beta,i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime})}}$

and

$\displaystyle
P^{\prime}_{(i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime},i\alpha
j\beta)}$ $\displaystyle=$ $\displaystyle\sum_{k\gamma
l\delta}P^{c}_{v,(i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime},k\gamma
l\delta)}P_{{\rm spg}{(k\gamma l\delta,i\alpha j\beta)}}$ (78)
$\displaystyle=$ $\displaystyle\frac{1}{N}\sum_{l}P_{{\rm
spg}{(i^{\prime}\alpha^{\prime}l\beta^{\prime},i\alpha j\beta)}}$
$\displaystyle=$ $\displaystyle\frac{1}{N}\sum_{l}P_{{\rm spg}{(i\alpha
j\beta,i^{\prime}\alpha^{\prime}l\beta^{\prime})}}.$

The final equality is derived from the symmetric property of $\bm{P}_{\rm
spg}$. Using the definition of $\bm{P}_{\rm spg}$, these elements are
expressed as

$\displaystyle P^{\prime}_{(i\alpha
j\beta,i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime})}$
$\displaystyle=\frac{1}{N|\mathcal{G}|}\sum_{l}\sum_{\hat{g}\in\mathcal{G}}\Gamma(\hat{g})_{(i\alpha,i^{\prime}\alpha^{\prime})}\Gamma({\hat{g}})_{(l\beta,j^{\prime}\beta^{\prime})}$
(79)

and

$\displaystyle
P^{\prime}_{(i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime},i\alpha
j\beta)}$
$\displaystyle=\frac{1}{N|\mathcal{G}|}\sum_{l}\sum_{\hat{g}\in\mathcal{G}}\Gamma(\hat{g})_{(i\alpha,i^{\prime}\alpha^{\prime})}\Gamma({\hat{g}})_{(j\beta,l\beta^{\prime})},$
(80)

respectively. The summations over index $l$ in Eqs. (B) and (B) must then be
equal, i.e.,

$\sum_{l}\Gamma({\hat{g}})_{(l\beta,j^{\prime}\beta^{\prime})}=\sum_{l}\Gamma({\hat{g}})_{(j\beta,l\beta^{\prime})},$
(81)

because $\Gamma({\hat{g}})_{(l\beta,j^{\prime}\beta^{\prime})}$ can be non-
zero only for single atom $l$ such that operation $\hat{g}$ moves atom $l$ to
atom $j^{\prime}$ and its value corresponds to the representation of $\hat{g}$
for the Cartesian components $(\beta,\beta^{\prime})$. When we denote these
summations as $\Gamma^{\prime}(\hat{g})_{(\beta,\beta^{\prime})}$, the
elements of $\bm{P}^{\prime}$ are written as

$\displaystyle P^{\prime}_{(i\alpha
j\beta,i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime})}=P^{\prime}_{(i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime},i\alpha
j\beta)}$
$\displaystyle=\frac{1}{N|\mathcal{G}|}\sum_{\hat{g}\in\mathcal{G}}\Gamma(\hat{g})_{(i\alpha,i^{\prime}\alpha^{\prime})}\Gamma^{\prime}(\hat{g})_{(\beta,\beta^{\prime})},$
(82)

which is equivalent to $\bm{P}_{v}^{c}\bm{P}_{\rm
spg}=(\bm{P}_{v}^{c}\bm{P}_{\rm spg})^{\top}$. Therefore,
$\bm{P}_{v}\bm{P}_{\rm spg}=(\bm{P}_{v}\bm{P}_{\rm spg})^{\top}$ can also be
derived, and it has been proved that $[\bm{P}_{\rm spg},\bm{P}_{v}]=0$.

(2) Proof of $[\bm{P}_{\rm spg},\bm{P}_{w}]=0$. We demonstrate that
$\bm{P}_{w}\bm{P}_{\rm spg}$ is symmetric. The projection matrix $\bm{P}_{w}$
is represented by $\bm{P}_{w}=\bm{B}_{w}\bm{B}_{w}^{\top}$. Therefore, the
elements of matrix $\bm{P}^{\prime\prime}=\bm{P}_{w}\bm{P}_{\rm spg}$ are
written as

$P^{\prime\prime}_{(i\alpha
j\beta,i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime})}=\frac{1}{2}\left[P_{{\rm
spg}{(i\alpha
j\beta,i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime})}}+P_{{\rm
spg}{(j\beta
i\alpha,i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime})}}\right]$ (83)

and

$P^{\prime\prime}_{(i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime},i\alpha
j\beta)}=\frac{1}{2}\left[P_{{\rm
spg}{(i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime},i\alpha
j\beta)}}+P_{{\rm
spg}{(j^{\prime}\beta^{\prime}i^{\prime}\alpha^{\prime},i\alpha
j\beta)}}\right].$ (84)

Then, the projection matrix $\bm{P}_{\rm spg}$ is orthogonal and has the
following symmetric properties of

$P_{{\rm spg}{(i\alpha
j\beta,i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime})}}=P_{{\rm
spg}{(i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime},i\alpha j\beta)}}$
(85)

and

$\displaystyle P_{{\rm spg}{(j\beta
i\alpha,i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime})}}$ $\displaystyle=$
$\displaystyle P_{{\rm
spg}{(i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime},j\beta i\alpha)}}$
(86) $\displaystyle=$ $\displaystyle P_{{\rm
spg}{(j^{\prime}\beta^{\prime}i^{\prime}\alpha^{\prime},i\alpha j\beta)}},$

where the final equality is derived from the fact that $\bm{P}_{\rm spg}$ is
expressed as the sum of the Kronecker products of representation
$\bm{\Gamma}(\hat{R})$. After applying these equations, we can conclude that
$\bm{P}^{\prime\prime}$ is symmetric as

$P^{\prime\prime}_{(i\alpha
j\beta,i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime})}=P^{\prime\prime}_{(i^{\prime}\alpha^{\prime}j^{\prime}\beta^{\prime},i\alpha
j\beta)},$ (87)

which proves that $[\bm{P}_{\rm spg},\bm{P}_{w}]=0$.

##  Appendix C Efficient implementation of force-constant estimation

When dividing the displacement-force dataset, we adopt the following procedure
to compute $\bm{X}^{\top}\bm{X}$ and $\bm{X}^{\top}\bm{y}$ more efficiently.
Introducing variables $z_{i\alpha,c}^{(s,{\rm FC2})}$ and
$z_{i\alpha,c}^{(s,{\rm FC3})}$ defined by

$\displaystyle z_{i\alpha,c}^{(s,{\rm FC2})}$ $\displaystyle=$
$\displaystyle-\sum_{j\beta}u_{j\beta}^{(s)}C_{i\alpha,j\beta,c}^{\rm(FC2)}$
$\displaystyle z_{i\alpha,c}^{(s,{\rm FC3})}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{j\beta}\sum_{k\gamma}u_{j\beta}^{(s)}u_{k\gamma}^{(s)}C_{i\alpha,j\beta,k\gamma,c}^{\rm(FC3)},$
(88)

the forces acting on atoms in supercell structure $s$ is rewritten by

$\displaystyle f_{i\alpha}^{(s)}$ $\displaystyle=$
$\displaystyle\sum_{b}w_{b}^{\rm(FC2)}\sum_{c}z_{i\alpha,c}^{(s,{\rm
FC2})}E_{c,b}^{\rm(FC2)}$ (89) $\displaystyle+$
$\displaystyle\sum_{b}w_{b}^{\rm(FC3)}\sum_{c}z_{i\alpha,c}^{(s,{\rm
FC3})}E_{c,b}^{\rm(FC3)}.$

In a matrix form, the forces are represented as

$\displaystyle\bm{f}^{(s)}$ $\displaystyle=$
$\displaystyle\bm{Z}^{(s)}\bm{E}\bm{w}$ (90) $\displaystyle=$
$\displaystyle\bm{X}^{(s)}\bm{w},$

where

$\displaystyle\bm{Z}^{(s)}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}\bm{Z}^{(s,\rm FC2)}&\bm{Z}^{(s,\rm FC3)}\\\
\end{bmatrix}$ (91) $\displaystyle=$
$\displaystyle\begin{bmatrix}\cdots&\vdots&\cdots&\vdots&\cdots\\\
\cdots&z_{i\alpha,c}^{(s,{\rm FC2})}&\cdots&z_{i\alpha,c}^{(s,{\rm
FC3})}&\cdots\\\ \cdots&\vdots&\cdots&\vdots&\cdots\\\ \end{bmatrix},$

and

$\bm{E}=\begin{bmatrix}\bm{E}^{(\rm FC2)}&\bm{0}\\\ \bm{0}&\bm{E}^{(\rm
FC3)}\\\ \end{bmatrix}.$ (92)

For the displacement-force dataset, matrix $\bm{Z}$ is expressed as

$\displaystyle\bm{Z}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}\bm{Z}^{(1)}\\\ \bm{Z}^{(2)}\\\ \bm{Z}^{(3)}\\\
\vdots\\\ \bm{Z}^{(N_{s})}\end{bmatrix}=\begin{bmatrix}\bm{Z}^{(1,\rm
FC2)}&\bm{Z}^{(1,\rm FC3)}\\\ \bm{Z}^{(2,\rm FC2)}&\bm{Z}^{(2,\rm FC3)}\\\
\bm{Z}^{(3,\rm FC2)}&\bm{Z}^{(3,\rm FC3)}\\\ \vdots&\vdots\\\
\bm{Z}^{(N_{s},\rm FC2)}&\bm{Z}^{(N_{s},\rm FC3)}\\\ \end{bmatrix}$ (93)
$\displaystyle=$ $\displaystyle\begin{bmatrix}\bm{Z}^{(\rm FC2)}&\bm{Z}^{(\rm
FC3)}\\\ \end{bmatrix}.$

Using these representations, $\bm{X}^{\top}\bm{X}$ and $\bm{X}^{\top}\bm{y}$
can be rewritten as

$\displaystyle\bm{X}^{\top}\bm{X}$ $\displaystyle=$
$\displaystyle\bm{E}^{\top}\bm{Z}^{\top}\bm{Z}\bm{E}=\begin{bmatrix}\bm{E}^{(\rm
FC2)\top}\bm{Z}^{(\rm FC2)\top}\bm{Z}^{(\rm FC2)}\bm{E}^{(\rm
FC2)}&\bm{E}^{(\rm FC2)\top}\bm{Z}^{(\rm FC2)\top}\bm{Z}^{(\rm
FC3)}\bm{E}^{(\rm FC3)}\\\ \bm{E}^{(\rm FC3)\top}\bm{Z}^{(\rm
FC3)\top}\bm{Z}^{(\rm FC2)}\bm{E}^{(\rm FC2)}&\bm{E}^{(\rm
FC3)\top}\bm{Z}^{(\rm FC3)\top}\bm{Z}^{(\rm FC3)}\bm{E}^{(\rm FC3)}\\\
\end{bmatrix},$ (94)

and

$\displaystyle\bm{X}^{\top}\bm{y}$ $\displaystyle=$
$\displaystyle\bm{E}^{\top}\bm{Z}^{\top}\bm{y}=\begin{bmatrix}\bm{E}^{(\rm
FC2)\top}\bm{Z}^{(\rm FC2)\top}\bm{y}\\\ \bm{E}^{(\rm FC3)\top}\bm{Z}^{(\rm
FC3)\top}\bm{y}\\\ \end{bmatrix},$ (95)

respectively. Finally, the matrix products with forms of $\bm{Z}^{\top}\bm{Z}$
and $\bm{Z}^{\top}\bm{y}$ are calculated efficiently by dividing $\bm{Z}$ and
$\bm{y}$ into subsets.

## References

  * (1) A. Seko and A. Togo, “Symfc,” https://github.com/symfc/symfc.
  * Dove (1993) M. T. Dove,  _Introduction to Lattice Dynamics_ (Cambridge University Press, 1993).
  * Wallace (1998) D. C. Wallace,  _Thermodynamics of crystals_ (Dover Publications, 1998).
  * Ziman (1960) J. M. Ziman,  _Electrons and phonons: the theory of transport phenomena in solids_ (Oxford University Press, 1960).
  * Parlinski  _et al._ (1997) K. Parlinski, Z. Q. Li, and Y. Kawazoe, Phys. Rev. Lett. 78, 4063 (1997).
  * Hellman  _et al._ (2013) O. Hellman, P. Steneteg, I. A. Abrikosov, and S. I. Simak, Phys. Rev. B 87, 104111 (2013).
  * Tadano and Tsuneyuki (2018) T. Tadano and S. Tsuneyuki, J. Phys. Soc. Jpn. 87, 041015 (2018).
  * Eriksson  _et al._ (2019) F. Eriksson, E. Fransson, and P. Erhart, Adv. Theory Simul. 2, 1800184 (2019).
  * Togo (2023) A. Togo, J. Phys. Soc. Jpn. 92, 012001 (2023).
  * Errea  _et al._ (2013) I. Errea, M. Calandra, and F. Mauri, Phys. Rev. Lett. 111, 177002 (2013).
  * van Roekeghem  _et al._ (2021) A. van Roekeghem, J. Carrete, and N. Mingo, Comput. Phys. Commun. 263, 107945 (2021).
  * Nye (1985) J. F. Nye,  _Physical Properties of Crystals_ (Oxford University Press, 1985).
  * El-Batanouny and Wooten (2008) M. El-Batanouny and F. Wooten,  _Symmetry and Condensed Matter Physics: A Computational Approach_ (Cambridge University Press, 2008).
  * Venkataraman  _et al._ (1975) G. Venkataraman, L. A. Feldkamp, and V. C. Sahni,  _Dynamics of perfect crystals_ (MIT press, 1975).
  * Zhou  _et al._ (2014) F. Zhou, W. Nielson, Y. Xia, and V. Ozoliņš, Phys. Rev. Lett. 113, 185501 (2014).
  * Tadano and Tsuneyuki (2015) T. Tadano and S. Tsuneyuki, Phys. Rev. B 92, 054301 (2015).
  * Yanai  _et al._ (2011) H. Yanai, K. Takeuchi, and Y. Takane,  _Projection Matrices, Generalized Inverse Matrices, and Singular Value Decomposition_, Statistics for social and behavioral sciences (Springer, Dordrecht, 2011).
  * Strang (2023) G. Strang,  _Introduction to Linear Algebra_ , 6th ed. (CUP, 2023).
  * Piziak  _et al._ (1999) R. Piziak, P. Odell, and R. Hahn, Comput. Math. Appl. 37, 67 (1999).
  * von Neumann (1949) J. von Neumann, Annals of Math. 50, 401 (1949).
  * Halperin (1962) I. Halperin, Acta Sci. Math.(Szeged) 23, 96 (1962).
  * Galántai (2013) A. Galántai,  _Projectors and projection methods_ , Vol. 6 (Springer Science & Business Media, 2013).
  * Sedgewick (2001) R. Sedgewick,  _Algorithms in c, part 5: graph algorithms, third edition_ , 3rd ed. (Addison-Wesley Professional, 2001).
  * Seko  _et al._ (2019) A. Seko, A. Togo, and I. Tanaka, Phys. Rev. B 99, 214108 (2019).
  * Seko (2023) A. Seko, J. Appl. Phys. 133, 011101 (2023).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Peierls (1929) R. E. Peierls, Ann. Phys. 395, 1055 (1929).
  * Peierls (2001) R. E. Peierls,  _Quantum theory of solids_ (Oxford University Press, 2001).
  * Allen and Perebeinos (2018) P. B. Allen and V. Perebeinos, Phys. Rev. B 98, 085427 (2018).
  * Togo  _et al._ (2015) A. Togo, L. Chaput, and I. Tanaka, Phys. Rev. B 91, 094306 (2015).
  * Togo  _et al._ (2023) A. Togo, L. Chaput, T. Tadano, and I. Tanaka, J. Phys. Condens. Matter 35, 353001 (2023).
  * Fujii and Seko (2022) S. Fujii and A. Seko, Comput. Mater. Sci. 204, 111137 (2022).
  * (37) A. Seko, Polynomial Machine Learning Potential Repository at Kyoto University, https://sekocha.github.io.
  * Togo and Seko (2024) A. Togo and A. Seko, arXiv preprint arXiv:2401.17531 (2024).