# Fileset

[2403.03588v3.pdf](https://mdr.nims.go.jp/filesets/acf6d486-498f-49cd-aa2e-1888c04314c6/download)

## Creator

[Atsuto Seko](https://orcid.org/0000-0002-2473-3837), [Atsushi Togo](https://orcid.org/0000-0001-8393-9766)

## Rights

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

## Other metadata

[Projector-based efficient estimation of force constants](https://mdr.nims.go.jp/datasets/a11b568c-34e2-4306-bcb5-169417d3b730)

## Fulltext

Projector-based efficient estimation of force constantsAtsuto Seko1, ∗ and Atsushi Togo21Department of Materials Science and Engineering, Kyoto University, Kyoto 606-8501, Japan2Center for Basic Research on Materials National Institute for Materials Science, Tsukuba, Ibaraki 305-0047, Japan(Dated: October 29, 2024)Estimating force constants for crystal structures is crucial for calculating various phonon-relatedproperties. However, this task becomes particularly challenging when dealing with a large number ofatoms or when third- and higher-order force constants are required. In this study, we propose an ef-ficient approach that involves constructing a complete orthonormal basis set for the force constants.We formulate this approach using projection matrices and their eigenvectors to meet the require-ments for the force constants. This basis set enables precise inference of the force constants fromdisplacement-force datasets. Our efficient algorithms for basis-set construction and force constantestimation are implemented in the symfc code. Furthermore, several applications demonstrated inthis study indicate that the current approach facilitates the efficient and accurate determination offorce constants.I. INTRODUCTIONForce constants in crystal structures are essential forcalculating various phonon-related properties [1–3]. Toestimate these force constants, supercell approaches com-bined with the finite displacement method are commonlyemployed [4–8]. Additionally, in self-consistent phononapproaches, effective force constants are often derivedfrom supercell structures that replicate the atomic dis-tribution at specific temperatures [9, 10]. These meth-ods involve introducing random or systematic atomic dis-placements into supercells, calculating the forces on theatoms, and then estimating the supercell force constantsfrom the resulting displacement-force dataset.Determining the force-constant elements requires nu-merous supercells with different displacement configura-tions, because the number of elements to be determinedcan be substantial. At the same time, the number offorce-constant elements can be reduced by using a basisset that represents the force constants, applying permu-tation symmetry rules, which stem from the definitionof force constants as derivatives of the potential energy,and sum rules, which derive from the invariance under in-finitesimal translations [2]. Additionally, the symmetricproperties of the structure can further reduce the num-ber of required force-constant elements [6, 7, 11]. A sig-nificantly reduced but complete orthonormal basis setenables efficient and accurate calculation of force con-stants that satisfy necessary requirements [12, 13]. How-ever, constructing such a basis set and estimating theforce constants using the basis set can be computation-ally intensive, particularly for systems with many atomsor when higher-order force constants are involved.Various approaches, including those utilizing com-pressed sensing, have been developed to estimate forceconstants with high precision [14, 15]. Despite the avail-ability of other approaches [5–7], force constant estima-∗ seko@cms.mtl.kyoto-u.ac.jption remains a computationally challenging task. In thisstudy, we propose a projector-based approach that effi-ciently constructs a complete orthonormal basis set forforce constants. We introduce several compression tech-niques that allow for a reduction in the size of the ba-sis set. These techniques involve compressing a pro-jection matrix using eigenvectors of another projectionmatrix onto a subspace defined by partial requirementsfor the force constants. The current procedures for effi-cient basis-set construction and force constant estimationfrom displacement-force datasets are implemented in thesymfc code [16]. We demonstrate various applications ofthis approach, showing that it enables the efficient andprecise determination of force constants.Section II introduces a straightforward projector-basedformulation for constructing basis sets for force con-stants. Section III presents various techniques for effi-ciently constructing basis sets using projection matrices,enabling the construction of basis sets for higher-orderforce constants. Section III I outlines the current proce-dure for constructing basis sets that exactly satisfy therequired constraints for force constants. Section IV de-scribes the procedure for estimating force constants froma displacement-force dataset using the constructed basisset. Section V demonstrates several applications of theproposed projector-based procedures for basis set con-struction and force constant estimation in diamond sil-icon and ionic systems. Finally, Sec. VI provides theconcluding remarks.II. PROJECTION MATRICES FOR FORCECONSTANTSA. Crystal potential and force constantsForce constants of a crystal are derived from the Tay-lor expansion of the potential energy V with respect toatomic displacements. In this study, we consider super-cell force constants, which are estimated from a set offinite displacements and forces acting on atoms in super-arXiv:2403.03588v3  [cond-mat.mtrl-sci]  26 Oct 2024mailto:seko@cms.mtl.kyoto-u.ac.jp2cells. For a supercell with N atoms, the potential energyis represented byV = Θ0 +∑iαΘiαuiα +12∑iα,jβΘiα,jβuiαujβ+13!∑iα,jβ,kγΘiα,jβ,kγuiαujβukγ + · · · , (1)where i, j, and k are the atomic indices within the super-cell, and α, β, and γ are the Cartesian components. Thesecond-order and third-order supercell force constants arerepresented by Θiα,jβ and Θiα,jβ,kγ , respectively. In ad-dition, we assume that the first-order force constants Θiαis equal to zero because atomic displacements are mea-sured from their equilibrium positions.B. Projection matrix for second-order forceconstantsIn the current formulation, the second-order force con-stants are represented as a column vector with 9N2 ele-ments, denoted asθ(FC2) = [Θ1x,1x,Θ1x,1y,Θ1x,1z,Θ1x,2x, · · · ]⊤, (2)where the components Θiα,jβ are listed in ascending or-der of the indices (i, α, j, β). We start by considering theforce constants that are invariant with the space groupoperations G for the target structure. Since we are dis-cussing supercell force constants, the lattice translationsubgroup in G is limited to those of the lattice points in-side the supercell. An operation ĝ ∈ G acts on the forceconstant vector asĝθ(FC2) = θ(FC2) [Γ(ĝ)⊗ Γ(ĝ)] , (3)where Γ(ĝ) denotes the matrix representation of opera-tion ĝ for the first-order force constants {Θiα}. We usethe group-theoretical projection operator method [12] toreduce the Kronecker product of representations to theidentity irreducible representation. The projection ma-trix for the identity irreducible representation is givenasPspg =1|G|∑ĝ∈GΓ(ĝ)⊗ Γ(ĝ), (4)where |G| denotes the order of group G.Then, we introduce projection matrices to enforce theconstraints associated with the sum rules and permuta-tion symmetry rules for the force constants. The sumrule for (iα, β) and the permutation symmetry rule for(iα, jβ) are expressed as∑jΘiα,jβ = 0 (5)andΘiα,jβ −Θjβ,iα = 0, (6)respectively. These rules impose linear constraints onthe force constants, which are mutually independent.In matrix form, these constraints are represented asC⊤θ(FC2) = 0. Therefore, a force-constant basis set thatsatisfies these linear constraints spans the kernel of C⊤,which can be expressed asKer(C⊤) = {θ(FC2) ∈ R9N2|C⊤θ(FC2) = 0}. (7)Since Ker(C⊤) is the complementary subspace spannedby C, the orthogonal projection matrix onto Ker(C⊤) isgiven by [17, 18]Pperm∩sum = I −C(C⊤C)−1C⊤, (8)where I denotes the unit matrix of size 9N2. A derivationof the projection matrix can be found in Appendix A.Each column ofC corresponds to either the sum rule orthe permutation symmetry rule for a specific combinationof indices. The sum rule for (i, α, β) can be expressedusing the column unit vector viαβ such that v⊤iαβθ(FC2) =0. The vector viαβ contains N non-zero elements, andthe set of vectors {viαβ} is orthonormal. Similarly, thepermutation symmetry rule for the combination (iα, jβ)is described by the column unit vector wiα,jβ such thatw⊤iα,jβθ(FC2) = 0. The vector wiα,jβ has only two non-zero elements: 1/√2 and −1/√2. Additionally, the set ofvectors {wiα,jβ} is orthonormal if the index combinationsare distinct and do not contain duplicates.Using a complete set of viαβ and wiα,jβ , the matrix Cis described asC = [Cv,Cw], (9)whereCv = [· · · ,viαβ , · · · ], Cw = [· · · ,wiα,jβ , · · · ]. (10)The matrices Cv and Cw have dimensions (9N2, 9N)and approximately (9N2, 9N2/2), respectively. Since thevectors viαβ and wiα,jβ are linearly independent butnot necessarily orthogonal, the projection matrix ontoKer(C⊤) is given by Eq. (8), which involves the inverseof C⊤C.Finally, we present the full projection matrix onto theintersection of vector subspaces determined by projectionmatrices Pspg and Pperm∩sum. Since projection matricesPspg and Pperm∩sum are commutative, as shown in Ap-pendix B, the full projection matrix can be written asP = PspgPperm∩sum [18, 19]. By solving the eigenvalueproblem for the full projection matrix:Pb = b, (11)we can obtain an orthonormal basis set for representingthe force constants, B. This basis set is invariant underspace group operations and satisfies both the sum rulesand the permutation symmetry rules. Since all non-zeroeigenvalues of the projection matrix are equal to one, thenumber of independent basis vectors corresponds to thetrace of the projection matrix.3C. Projection matrix for third-order forceconstantsThe third-order force constants are represented in avector form with 27N3 elements asθ(FC3) = [Θ1x,1x,1x,Θ1x,1x,1y,Θ1x,1x,1z,Θ1x,1x,2x, · · · ]⊤,(12)which are listed in ascending order of (i, α, j, β, k, γ). Anoperation ĝ ∈ G acts on the third-order force constantvector asĝθ(FC3) = θ(FC3) [Γ(ĝ)⊗ Γ(ĝ)⊗ Γ(ĝ)] . (13)Therefore, the projection matrix for the identity irre-ducible representation is given asPspg =1|G|∑ĝ∈GΓ(ĝ)⊗ Γ(ĝ)⊗ Γ(ĝ). (14)The sum rules and permutation symmetry rules can alsobe formulated in matrix form as C⊤θ(FC3) = 0, analo-gous to the matrix representation used for second-orderforce constants. The kernel of C⊤ is written asKer(C⊤) = {θ(FC3) ∈ R27N3|C⊤θ(FC3) = 0}, (15)and the orthogonal projection matrix onto Ker(C⊤) isgiven by [17, 18]Pperm∩sum = I −C(C⊤C)−1C⊤, (16)where I denotes the unit matrix of size 27N3.The sum rule for (i, α, j, β, γ) is given by∑kΘiα,jβ,kγ = 0, (17)and can be expressed as v⊤iαjβγθ(FC3) = 0. The vectorviαjβγ contains only N non-zero elements out of a to-tal of 27N3 elements, and the set of vectors {viαjβγ} ismutually orthonormal.Representing the permutation symmetry rules forthird-order force constants is more complex than forsecond-order force constants. For the permutation sym-metry of (iα, jβ, kγ), the rules are expressed asΘiα,jβ,kγ −Θiα,kγ,jβ = 0Θiα,jβ,kγ −Θjβ,iα,kγ = 0Θiα,jβ,kγ −Θjβ,kγ,iα = 0 (18)Θiα,jβ,kγ −Θkγ,iα,jβ = 0Θiα,jβ,kγ −Θkγ,jβ,iα = 0,when all indices iα, jβ, and kγ are distinct. Theseconstraints can be represented using a matrix com-posed of column unit vectors Wiα,jβ,kγ such thatW⊤iα,jβ,kγθ(FC3) = 0. Each column vector of Wiα,jβ,kγcontains two non-zero elements, 1/√2 and −1/√2. Al-though these vectors are not orthogonal, they are linearlyindependent. For the case where the indices are of theform (iα, iα, kγ), the permutation symmetry rules aregiven byΘiα,iα,kγ −Θiα,kγ,iα = 0Θiα,iα,kγ −Θkγ,iα,iα = 0. (19)Similarly, the matrix of column unit vectors Wiα,iα,kγis defined such that W⊤iα,iα,kγθ(FC3) = 0. Each columnvector in Wiα,iα,kγ has two non-zero elements, 1/√2 and−1/√2.Using a linearly independent complete set of viαjβγand Wiα,jβ,kγ , the matrix C is described asC = [Cv,Cw], (20)whereCv = [· · · ,viαjβγ , · · · ], Cw = [· · · ,Wiα,jβ,kγ , · · · ]. (21)The matrices Cv and Cw have dimensions (27N3, 27N2)and approximately (27N3, (5/6)27N3), respectively. Fi-nally, the full projection matrix is given by P =PspgPperm∩sum. By solving the eigenvalue problem forthe full projection matrix, Pb = b, an orthonormal ba-sis set for representing third-order force constants is ob-tained. Note that the derivation of basis sets for rep-resenting higher-order force constants can be similarlyformalized. The symfc code implements the constructionof basis sets up to fourth-order force constants.III. EFFICIENT BASIS SET CONSTRUCTIONThe previous section outlines a straightforward formu-lation of the projector-based approach for constructing acomplete basis set for force constants. However, a directapplication of this method requires computing the inverseof C⊤C and solving large eigenvalue problems. Specif-ically, for second-order force constants, the eigenvalueproblem is of size (9N2, 9N2), and for third-order forceconstants, it is of size (27N3, 27N3). Such computationscan be prohibitively expensive, particularly for the third-order force constants. In Sec. III A–IIIH, we introduceand formulate various techniques designed to efficientlyconstruct a basis set for force constants, thereby address-ing associated computational challenges. The overall pro-cedure for constructing a basis set is outlined in Sec. III I.A. Decomposition of projection matrixInstead of deriving the full projection matrix by com-puting the inverse ofC⊤C, we construct a complete basisset using an approximation of the full projection matrix.This approximated projection matrix is composed of in-dividual projection matrices that apply the sum rules andpermutation symmetry rules separately. For the second-order force constants, orthogonal projection matrices Pv4and Pw onto Ker(C⊤v ) and Ker(C⊤w ) are expressed usingorthonormal vector sets {viαβ} and {wiα,jβ} as [17]Pv = I −N∑i=13∑α=13∑β=1viαβv⊤iαβ= I −CvC⊤v (22)andPw = I −∑(iα,jβ)wiα,jβw⊤iα,jβ= I −CwC⊤w , (23)respectively. For the third-order force constants, orthog-onal projection matrices Pv and Pw onto Ker(C⊤v ) andKer(C⊤w ) are given asPv = I −N∑i=13∑α=1N∑j=13∑β=13∑γ=1viαjβγv⊤iαjβγ= I −CvC⊤v (24)andPw = I −Cw(C⊤wCw)−1C⊤w , (25)respectively.Thus, projection matrices that apply the sum rulesand permutation symmetry rules can be defined indi-vidually. However, the projection matrix Pperm∩sum isnot simply the product of Pw and Pv due to the non-commutative nature of these matrices. The projectionmatrix Pperm∩sum along with its eigenvectors will be de-rived as demonstrated in Sec. III B.B. Non-commutative projection matricesGiven non-commutative projection matrices P1 andP2, the projection matrix P1∩2 can be expressed usingthe von Neumann–Halperin theorem as [19–22]P1∩2 = limp→∞P1(P2P1)p, (26)which represents the projection matrix onto the intersec-tion of the subspaces associated with P1 and P2. Thegeometric interpretation of Eq. (26) is that a vector xis first projected onto the subspace associated with P1,and the resulting vector is subsequently projected ontothe subspace associated with P2. This process of alter-nating projections onto the subspaces associated with P1and P2 is repeated. In the end, the projected vector con-verges to P1∩2x.For practical purposes, P1∩2 can be approximated us-ing a finite value of p. In this study, setting p = 1 and ap-proximating P1∩2 as P̃1∩2 ≃ P1(P2P1)1 suffice to derivea complete basis set for P1∩2. By solving the eigenvalueproblem for the approximated projection matrix definedasP̃1∩2 = P1P2P1 (27)and removing the eigenvectors with eigenvalues less thanone, a complete basis set for P1∩2 can be obtained. Thiscan be derived from the fact that P1∩2 = limp→∞ P̃ p1∩2.Note that the matrix described by Eq. (27) is an ap-proximation and does not fully satisfy the conditions ofa projection matrix. Nonetheless, the basis set formed byits eigenvectors corresponding to eigenvalues of one con-stitutes a complete set of eigenvectors for the projectionmatrix P1∩2.For the non-commutative matrices Pv and Pw, the pro-jection matrix Pperm∩sum and the full projection matrixP = PspgPperm∩sum can be approximated asP̃perm∩sum = PwPvPw (28)andP̃ = PspgPwPvPw, (29)respectively. Solving the eigenvalue problem for theapproximated full projection matrix and removing theeigenvectors with eigenvalues less than one yields a com-plete basis set for the full projection matrix.C. Compression of projection matricesThe compression of projection matrices is crucial forconstructing a complete force-constant basis set of thefull projection matrix without relying on approximations.Although the complete basis set can be derived by solv-ing the eigenvalue problem associated with the product ofpartial projection matrices given by Eq. (29), the compu-tational demands of full-sized large projection matricesfor third-order force constants, scaling as (27N3, 27N3),and higher-order force constants present significant chal-lenges for several reasons: (1) The construction of full-sized partial projection matrices is computationally in-tensive and necessitates substantial memory allocation.(2) The computation of products of full-sized partial pro-jection matrices is also computationally demanding. (3)Solving eigenvalue problems for full-sized projection ma-trices is prohibitively expensive in terms of computa-tional resources.Therefore, we propose a technique to compress projec-tion matrices, enabling the execution of these operationswith the smallest possible compressed projection matri-ces. This compression procedure is designed to retain allforce-constant components essential for constructing thecomplete set of the full projection matrix. It relies on theeigenvalue decompositions of partial projection matrices,where the eigenvectors can significantly enhance the com-pression of large matrices, provided that the eigenvec-tors of the partial projection matrix are well-compressed,5sparse, and efficiently computed. In the following, we willdemonstrate how to compress a projection matrix usingthe eigenvectors of a partial projection matrix.Firstly, we begin with the eigenvalue decompositionof a partial projection matrix. In this decomposition, apartial projection matrix Pj is expressed in terms of itscomplete and orthonormal set of eigenvectors Bj asPj = BjB⊤j , (30)whereBj has dimensions (9N2, nj) for second-order forceconstants and (27N3, nj) for third-order force constants.The variable nj denotes the number of columns, corre-sponding to the size of the complete set of eigenvectorsfor Pj . Since the eigenvectors Bj are well-compressed forthe partial projection matrices considered in this study,meaning that the size of the complete set of eigenvectorsis significantly smaller than the number of rows, they canbe utilized for compressing other projection matrices.The matrix Pi can be effectively compressed using theeigenvectors of Pj , denoted as Bj , when the productPjPiPj can be represented as a component of the fullprojection matrix. This compression is applicable whenthe eigenvectors are either known or readily computable.The compression process can be described as follows:PjPiPj = BjB⊤j PiBjB⊤j (Pj = BjB⊤j )= BjPi↓jB⊤j , (31)where Pi↓j represents the compressed version of Pi in thebasis defined by Bj . It is defined asPi↓j = B⊤j PiBj . (32)The size of Pi↓j is (nj , nj), which is significantly smallerthan the size of Pi. It is important to note that a com-plete set of eigenvectors Bj is employed to derive Pi↓j ,ensuring that this compression is exact and does not omitany basis vectors relevant to PjPiPj .The eigenvectors of PjPiPj can be efficiently com-puted by solving the eigenvalue problem for the smallermatrix Pi↓j as follows. (1) If the eigenvectors of Pi =BiB⊤i are known without explicitly constructing Pi andwithout solving its eigenvalue problem, then Pi↓j can bewritten asPi↓j = B⊤j BiB⊤i Bj . (33)Since B⊤j Bi is generally smaller in size, we first com-pute A = B⊤j Bi, and then find Pi↓j by calculatingPi↓j = AA⊤. If Pi must be constructed, Pi↓j can becomputed using Eq. (32). (2) The eigenvalue problemfor the compressed projection matrix Pi↓j is solved asPi↓j = Bi↓jB⊤i↓j , where Bi↓j are the compressed eigen-vectors by Bj . (3) The eigenvectors of PjPiPj witheigenvalues of one, Bi∩j , are determined asBi∩j = BjBi↓j , (34)because PjPiPj can be represented byPjPiPj = Bi∩jB⊤i∩j= BjBi↓jB⊤i↓jB⊤j . (35)Thus, a complete set of eigenvectors of PjPiPj can becalculated efficiently without the need to directly com-pute the product PjPiPj .If the projection matrices Pi and Pj commute, i.e.,[Pi,Pj ] = 0, this compression technique can be appliedto the product PiPj . Specifically, the product PiPj canbe expressed asPiPj = PiP2j = PjPiPj , (36)where we use the fact that projection matrices are idem-potent, P 2j = Pj .D. Application of sum rulesWe utilize the equation Pv = I − CvC⊤v , which ispresented in Eqs. (22) and (24), to incorporate the sumrules. The use of this formulation is computationally effi-cient because both the identity matrix I and the matrixCv are sparse. On the other hand, the eigenvectors ofPv, denoted by Bv, typically contain a larger number ofnon-zero elements compared to Cv.E. Application of permutation symmetry rulesThe eigenvectors of the projection matrix Pw, denotedby Bw, can be efficiently determined without the need toexplicitly construct Pw. These eigenvectors span the vec-tor space that satisfies the permutation symmetry rules.For second-order force constants, the permutation sym-metry rule for (iα, jβ) is specified by Eq. (23) and is rep-resented by a single vector wiα,jβ such that w⊤iα,jβθ = 0.Since the set of vectors {wiα,jβ} is orthonormal, the unitvector w⊥iα,jβ , which is orthogonal to wiα,jβ , can serveas an eigenvector of Pw. This vector w⊥iα,jβ has only twonon-zero elements of 1/√2 corresponding to the columns(i, α, j, β) and (j, β, i, α). In cases where iα = jβ, nopermutation symmetry rules are applicable. Therefore,w⊥iα,iα is defined as the vector with a single non-zero el-ement of one for the column (i, α, i, α). Thus, the eigen-vectors Bw are given byBw = [· · · ,w⊥iα,jβ , · · · ], (37)where the vectors {w⊥iα,jβ} are mutually orthonormal.The matrix Bw has dimensions (9N2, nw), where nw ≃9N2/2. Despite having 9N2 non-zero elements, the totalnumber of elements in Bw is much larger, highlightingthe sparsity of the eigenvectors.For third-order force constants, the permutation sym-metry rule for (iα, jβ, kγ) is specified by Eq. (18), which6comprises five distinct equations. This symmetry rulecan be represented in matrix form as W⊤iα,jβ,kγθ = 0,and the submatrix of Wiα,jβ,kγ , constructed by extract-ing only the non-zero rows, is given byWiα,jβ,kγ =1√21 1 1 1 1−1 0 0 0 00 −1 0 0 00 0 −1 0 00 0 0 −1 00 0 0 0 −1 . (38)Therefore, the column vector w⊥iα,jβ,kγ , which is or-thonormal to all column vectors of Wiα,jβ,kγ , can beconsidered as an eigenvector of Pw. It is given byw⊥iα,jβ,kγ =1√6[1 1 1 1 1 1]⊤, (39)which is represented in the submatrix form by omittingzero elements. When iα = jβ and iα ̸= kγ, the permuta-tion symmetry rule for (iα, iα, kγ) is given by Eq. (19).The two equations in Eq. (19) are also represented inmatrix form as W⊤iα,iα,kγθ = 0, where the submatrix ofWiα,iα,kγ comprising only non-zero rows is expressed asWiα,iα,kγ =1√2 1 1−1 00 −1 . (40)Consequently, the vector orthonormal to all column vec-tors of Wiα,jβ,kγ , denoted by w⊥iα,iα,kγ , is given byw⊥iα,iα,kγ =1√3[1 1 1]⊤. (41)If iα = jβ and iα = kγ, there are no permutation symme-try rules applicable. In this case, the vector w⊥iα,iα,iα isdefined as the vector with a single non-zero element equalto one for the column corresponding to (i, α, i, α, i, α).Thus, the eigenvectors Bw are written asBw = [· · · ,w⊥iα,jβ,kγ , · · · ], (42)where the vectors w⊥iα,jβ,kγ are mutually orthonormal.The matrix has dimensions (27N3, nw), where nw ≃27N3/6. Although the number of non-zero elements inBw is 27N3, it is significantly smaller than the total num-ber of elements in Bw.F. Effective use of lattice translation groupHere, we introduce the projection matrix correspond-ing to the identity irreducible representation of the lat-tice translation group, which is a normal subgroup of thespace group. While the inclusion of the projection matrixfor the lattice translation group is not strictly necessaryfor the comprehensive description of the full projectionmatrix, it is advantageous for the compression of the pro-jection matrices that constitute components of the fullprojection matrix. The projection matrix for the latticetranslation group can be derived using a methodologyanalogous to that employed for the projection matrix as-sociated with permutation symmetry rules.The force constants are invariant under operationsfrom the lattice translation group T = {t1, t2, · · · }. Here,representative atoms with respect to lattice translationsare denoted by indices {i′}, and the translation of atomi is represented as tn(i). For the second-order force con-stants involving a representative atom i′, the invarianceunder a lattice translation tn is expressed byΘi′α,jβ −Θtn(i′)α,tn(j)β = 0, (43)where tn is a non-identity translation operation. Con-sequently, for each pair (i′, α, j, β), there are |T | − 1 in-variance equations, which can be represented in matrixform by Ti′α,jβ such that T⊤i′α,jβθ = 0. The submatrixof Ti′α,jβ composed only of non-zero rows is expressed asTi′,α,jβ =1√21 1 1 1−1 0 0 00 −1 0 0 · · ·0 0 −1 00 0 0 −1.... . . , (44)where the size of this submatrix is (|T |, |T | − 1). Thevector orthonormal to all column vectors of Ti′α,jβ canbe an eigenvector of PLT. This eigenvector is representedin the submatrix form ast⊥i′α,jβ =1√|T |[1 1 1 · · · ]⊤. (45)The eigenvectors of PLT, denoted by BLT, are then givenbyBLT = [· · · , t⊥i′α,jβ , · · · ], (46)where the set of vectors {t⊥i′α,jβ} is orthonormal. The sizeof BLT is (9N2, 9N2/|T |), and it contains 9N2 non-zeroelements.The third-order force constants for triplet atoms in-cluding a representative atom i′ are also invariant underthe operations of the lattice translation group T . Thisinvariance is expressed asΘi′α,jβ,kγ −Θtn(i′)α,tn(j)β,tn(k)γ = 0, (47)where tn is a translation operation that is not the iden-tity. Similarly to the case of second-order force con-stants, there are |T | − 1 equations for each triplet ofindices (i′, α, j, β, k, γ), which can be represented in ma-trix form as T⊤i′α,jβ,kγθ(FC3) = 0. The submatrix formsof Ti′α,jβ,kγ and t⊥i′α,jβ,kγ are obtained by extending theindices from (i′α, jβ) to (i′α, jβ, kγ) in Eqs. (44) and7(45). Therefore, the eigenvectors of PLT for the third-order force constants can be expressed asBLT = [· · · , t⊥i′α,jβ,kγ , · · · ], (48)where the set of vectors {t⊥i′α,jβ,kγ} is orthonormal. Thesize of BLT is (27N3, 27N3/|T |), and the number of non-zero elements is 27N3, which is considerably smaller thanthe total number of elements in BLT. Thus, BLT remainsa sparse matrix.The projection matrix PLT for the lattice translationgroup is given by its eigenvectors as PLT = BLTB⊤LT.The projection matrix for the lattice translation groupsatisfies the following relationship with the projectionmatrix Pspg for the space group:PLTPspg = PspgPLT = Pspg. (49)This equation holds because we have tnG = G for anytranslational operation tn ∈ T that is also included inthe group G. Therefore, the projection matrix for thelattice translation group leaves Pspg invariant.G. Properties of projection matricesThe commutation relations between different projec-tion matrices play a crucial role in the efficient construc-tion of force-constant basis sets. These relations include[Pspg,Pv] = 0, [Pspg,Pw] = 0, [Pv,Pw] ̸= 0. (50)Additionally, the projection matrices for the lattice trans-lation group PLT exhibit useful commutation propertiesexpressed as[PLT,Pv] = 0, [PLT,Pw] = 0, (51)in conjunction with the relationship of Eq. (49). In ad-dition, all these projection matrices are idempotent andsymmetric, which are represented byP 2spg = Pspg, P2v = Pv, P2w = Pw, P2LT = PLT (52)andP⊤spg = Pspg, P⊤v = Pv, P⊤w = Pw, P⊤LT = PLT, (53)respectively.H. Eigenvalue problems for projection matricesSolving eigenvalue problems for large and sparse pro-jection matrices is a fundamental element of the cur-rent procedure. As demonstrated in Sec. III I, the ef-ficiency of the solver is crucial in determining the designof the procedure. When addressing eigenvalue problemsfor sparse projection matrices, exploiting the underly-ing block-diagonal structures can be particularly advan-tageous.We consider an orthogonal projection matrix expressedin block-diagonal form, which is transformed by row andcolumn permutations. The block-diagonal form of theorthogonal projection matrix is given byM⊤PM =I 0 00 J 00 0 K , (54)where M denotes the permutation matrix for columns.The submatrices I, J , and K are all orthogonal pro-jection matrices, ensuring that P is also an orthogonalprojection matrix. When solving eigenvalue problems forthese submatrices, and obtaining the eigenvectors of I,J , and K asIbI = bI , JbJ = bJ , KbK = bK , (55)the vectorsb =bI00 , 0bJ0 , 00bK (56)can be adopted as eigenvectors of M⊤PM . The setof eigenvectors {b} is orthonormal, and the trace ofM⊤PM is given bytr(M⊤PM) = tr(I) + tr(J) + tr(K), (57)which equals the number of eigenvectors of M⊤PM .Consequently, the set of eigenvectors {b} constructedfrom {bI}, {bJ}, and {bK} forms a complete basis. Fi-nally, the eigenvectors of P can be computed as Mb.Thus, the eigenvectors of the original projection matrixP can be determined by solving the eigenvalue problemsfor the smaller matrices.We use a graph-based approach to determine the per-mutation matrix M and reveal the block-diagonal struc-ture of P . We represent the orthogonal projection matrixP as an undirected graph, where each row and columncorresponds to a vertex. Vertices s and t are adjacentif the element in P at row s and column t is non-zero.We then apply an efficient graph algorithm to identifythe connected components [23] of this graph. The algo-rithm groups the rows and columns of P into connectedcomponents, which correspond to block structures in thematrix. Finally, we construct the permutation matrix Mso that columns within the same connected componentare placed adjacent to each other in the reordered matrixM⊤PM .I. Procedure to calculate eigenvectors of fullprojection matrixAs described in Sec. III B, the eigenvectors of the ap-proximated projection matrix P̃ = PspgPwPvPw corre-sponding to eigenvalues equal to one constitute a com-plete basis for the full projection matrix P . The compu-tational steps for determining the eigenvectors of P areoutlined as follows.8(1) The projection matrix Pw is compressed usingBLT.The compressed projection matrix is given byPw↓LT = B⊤LTBwB⊤wBLT. (58)Thus, the product PspgPw can be rewritten asPspgPw = PspgPLTPwPLT= PspgBLTB⊤LTBwB⊤wBLTB⊤LT (59)= PspgBLTPw↓LTB⊤LT,where the properties of PLT and Eq. (33) are utilized inthe derivation.(2) The eigenvalue problem for Pw↓LT is solved. Thecompressed projection matrix Pw↓LT is decomposed asPw↓LT = Bw↓LTB⊤w↓LT. (60)Using this decomposition, the product of projection ma-trices PLTPwPLT is decomposed asPLTPwPLT = BLTBw↓LTB⊤w↓LTB⊤LT= Bw∩LTB⊤w∩LT, (61)whereBw∩LT = BLTBw↓LT. (62)The set of vectors Bw∩LT spans the intersection of thesubspaces spanned by Bw and BLT.(3) The matrix Pspg is compressed using Bw∩LT asfollows:Pspg↓(w∩LT) = B⊤w∩LTPspgBw∩LT. (63)Using this compressed projection matrix, the productPspgPw is expressed asPspgPw = PLTPwPLTPspgPLTPwPLT= Bw∩LTB⊤w∩LTPspgBw∩LTB⊤w∩LT (64)= Bw∩LTPspg↓(w∩LT)B⊤w∩LT.(4) The eigenvalue problem for Pspg↓(w∩LT) is solved,and it is decomposed asPspg↓(w∩LT) = Bspg↓(w∩LT)B⊤spg↓(w∩LT). (65)Consequently, the product PspgPw is expressed asPspgPw = Bw∩LTBspg↓(w∩LT)B⊤spg↓(w∩LT)B⊤w∩LT= Bspg∩wB⊤spg∩w, (66)whereBspg∩w = Bw∩LTBspg↓(w∩LT). (67)(5) The projection matrix Pv is compressed usingBspg∩w asPv↓(spg∩w) = B⊤spg∩wPvBspg∩w. (68)Using this compressed form of Pv, the approximated pro-jection matrix P̃ is represented byP̃ = PspgPwPvPw= PspgPwPvPspgPw= Bspg∩wB⊤spg∩wPvBspg∩wB⊤spg∩w= Bspg∩wPv↓(spg∩w)B⊤spg∩w. (69)(6) The eigenvalue problem for Pv↓(spg∩w) is solved. Byeliminating the eigenvectors associated with eigenvaluesless than one, the full projection matrix P is obtained asP = Bspg∩wBv↓(spg∩w)B⊤v↓(spg∩w)B⊤spg∩w= BB⊤, (70)where the basis set for the force constants is given byB = Bspg∩wBv↓(spg∩w). (71)J. Practical approximationsThe computational bottleneck in the current procedurearises from solving the eigenvalue problem specified byEq. (70), particularly as the size of the basis set in-creases. Empirical observations indicate that the eigen-value problem becomes computationally intensive whenthe dimension of the square projection matrix Pv↓(spg∩w)exceeds 20000. To alleviate this issue, an approximationcan be employed to reduce the basis set by assigning zerovalues to a subset of force constant elements that are ex-pected to be near zero. Specifically, when third-orderforce constants with indices S are set to zero, the projec-tion matrix designed to eliminate these zero elements isgiven byPNZ = I −∑(iα,jβ,kγ)∈Sηiα,jβ,kγη⊤iα,jβ,kγ , (72)where ηiα,jβ,kγ represents a vector with a single non-zeroelement equal to one, corresponding to each element inS. Consequently, the projection matrix PNZ is a diagonalmatrix with zero entries for indices in S and one entriesfor all other indices. Since PNZ and a projection matrixP1 are generally non-commutative, i.e., [PNZ,P1] ̸= 0,a basis set spanning the intersection of vector subspacesdetermined by the projection matrices can be obtainedby solving the eigenvalue problem for PNZP1PNZ anddiscarding eigenvectors with eigenvalues less than one,as shown in Sec. III B.In practical terms, force constant elements are set tozero by eliminating the basis vectors corresponding toS from Bw and BLT. Specifically, when Θiα,jβ,kγ isset to zero, the column basis vectors in Bw and BLTthat contain non-zero elements associated with the row(iα, jβ, kγ) are removed. This means that projection ma-trices are applied to Pw and PLT as follows:Pw∩NZ =[I −∑ibw,ib⊤w,i]Pw (73)9andPLT∩NZ =[I −∑ibLT,ib⊤LT,i]PLT, (74)where bw,i and bLT,i represent the basis vectors removedfrom Bw and BLT, respectively. These projection ma-trices can be utilized by substituting Pw and PLT in theformulation of basis set construction.A systematic approach to introducing zero elementsinvolves employing a cutoff distance. The force constantΘiα,jβ,kγ is set to zero if the distance between any pairof atoms i, j, and k exceeds the cutoff distance. In thisstudy, the distance between two atoms within a givenperiodic supercell is defined as the minimum distanceamong those calculated using both the atoms in the su-percell and the replica atoms located in adjacent super-cells.IV. FORCE CONSTANT ESTIMATIONA. Force-constant basis set expansionThe forces acting on atoms {f (s)iα } in a supercell struc-ture s corresponding to atomic displacements {u(s)iα } iswritten asf(s)iα = −∑jβΘiα,jβu(s)jβ −12∑jβ∑kγΘiα,jβ,kγu(s)jβ u(s)kγ +· · · .(75)The force constants are then expanded in basis sets givenby Eq. (71) asΘiα,jβ =∑bc(FC2)b B(FC2)iα,jβ,b (76)andΘiα,jβ,kγ =∑bc(FC3)b B(FC3)iα,jβ,kγ,b, (77)where cb denotes the expansion coefficient for b-th basis.Using these basis expansions, the forces in structure s areexpressed asf(s)iα = −∑bc(FC2)b∑jβu(s)jβ B(FC2)iα,jβ,b− 12∑bc(FC3)b∑jβ∑kγu(s)jβ u(s)kγB(FC3)iα,jβ,kγ,b+ · · · , (78)and they are rewritten asf(s)iα =∑bc(FC2)b x(s,FC2)iα,b +∑bc(FC3)b x(s,FC3)iα,b + · · · ,(79)wherex(s,FC2)iα,b = −∑jβu(s)jβ B(FC2)iα,jβ,b (80)andx(s,FC3)iα,b = −12∑jβ∑kγu(s)jβ u(s)kγB(FC3)iα,jβ,kγ,b. (81)Thus, when modeling forces up to third-order force con-stants, the forces in structure s are represented in matrixform as given byf (s) = X(s)c, (82)wheref (s) =...f(s)iα... , (83)X(s) =[X(s,FC2) X(s,FC3)]=· · ·... · · ·... · · ·· · · x(s,FC2)iα,b · · · x(s,FC3)iα,b · · ·· · ·... · · ·... · · · , (84)andc =[c(FC2)c(FC3)]=...c(FC2)b...c(FC3)b.... (85)B. Force constant estimationThe expansion coefficients c are determined from theforces acting on atoms in a training displacement-forcedataset consisting of Ns structures using standard linearregression. The training dataset is expressed in matrixform asX =X(1)X(2)X(3)...X(Ns) , y =f (1)f (2)f (3)...f (Ns) , (86)and the normal equations to find the least squares solu-tion are given byX⊤Xc = X⊤y. (87)Once the coefficients c are estimated, the force constantscan then be recovered from these coefficients usingθ(FC2) = B(FC2)c(FC2)θ(FC3) = B(FC3)c(FC3). (88)10Here, we consider solving the normal equations di-rectly. However, computing the entire matrix X can beresource-intensive, particularly with large displacement-force datasets and extensive basis sets. To optimize mem-ory usage and efficiently estimate the expansion coeffi-cients, we split the computations of X⊤X and executethem sequentially. This approach is based on the recur-sive method for linear regression [24]. When the datasetis partitioned into nq subsets, X⊤X and X⊤y can becalculated asX⊤X =nq∑q=1X⊤q Xq, X⊤y =nq∑q=1X⊤q yq, (89)where Xq and yq represent the q-th subsets of X and y,respectively. This approach eliminates the need to storethe entire matrix X in memory. The coefficients are thendetermined by solving the normal equations using fastlinear algebra libraries. A more detailed implementationfor efficiently computing X⊤X and X⊤y is provided inAppendix C.Note that the methodology for estimating higher-orderforce constants can be similarly formalized. The symfccode is capable of performing force constant estimationup to the fourth order.V. RESULTS AND DISCUSSIONA. Basis set construction1. Complete basis setFor various structures and supercell sizes, orthonor-mal basis sets for the representation of second-order andthird-order supercell force constants are determined us-ing the current projector-based approach implemented inthe symfc code. Given a unit cell structure and a spec-ified supercell size, a basis set for the force constants isderived by applying the procedure shown in Sec. III I.This basis set ensures that the supercell force constantscompletely satisfy the permutational symmetry require-ments, sum rules, and the symmetric properties of thesupercell.Figure 1 (a) illustrates the size of the basis set for thesecond-order force constants of various crystal structuresfor elemental systems. The total number of second-orderforce constant elements is 9N2, but this number is sig-nificantly reduced by applying constraints on the forceconstants. For second-order force constants, it is possibleto derive complete basis sets for supercells containing upto 10000 atoms and approximately 109 force constant el-ements, while still maintaining acceptable computationalefficiency.The total number of third-order force-constant ele-ments is 27N3, which can become quite substantial whenevaluating force constants for supercells of typical size.Consequently, the reduction of third-order force con-stants is significantly more computationally demandingcompared to that of second-order force constants for agiven N . The current projector-based approach, how-ever, efficiently reduces the number of force constantsand constructs orthonormal basis sets.Figures 1 (b) and (c) display the number of basis vec-tors required for third-order force constants in elementaland ionic structures, respectively. In these figures, thecutoff distance is not applied, and no force-constant el-ements are constrained to be zero. The computationalbottleneck in the current procedure arises from solvingthe eigenvalue problem given by Eq. (70), particularlywhen the basis set size becomes large. Thanks to thesmall primitive cells and high symmetry of elementalstructures, the current procedure can be applied to super-cells with up to 900 atoms using a standard workstation.In contrast, for ionic structures, the number of atomsin applicable supercells is reduced due to decreased lat-tice translations. For the well-known structures shownin Fig. 1 (c), the procedure generally performs well forsupercells containing up to 300 atoms. Additional detailson the computational costs associated with constructingcomplete basis sets are provided in Appendix D.2. Application of cutoff distanceWhen dealing with a large number of force constantelements or basis sets, it is practical to introduce a cutoffdistance to assign zero values to force constant elementsbeyond this distance. Figure 2 illustrates the dependenceof the number of basis vectors on the cutoff distance forthe diamond-type structure. The implementation of acutoff distance significantly reduces the number of ba-sis vectors, thereby alleviating the computational burdenassociated with basis set construction and force constantestimation.B. Force constant estimationWe will present various applications of the current ap-proach implemented in the symfc code for estimatingforce constants. The first set of applications involvescalculating lattice thermal conductivities (LTCs) fromdisplacement-force datasets generated using density func-tional theory (DFT) calculations. The second set focuseson calculating LTCs using machine learning potentials(MLPs). These applications demonstrate that the forceconstants can be estimated efficiently with the currentapproach. Furthermore, they emphasize the importanceof selecting appropriate parameters and provide recom-mended settings for estimating both second-order andthird-order force constants using this approach.111001021041061081010 0  1000  2000  3000  4000  5000Number of basis vectorsNumber of atomsBCCFCCHCPSCDiamond1001021041061081010 0  300  600  900Number of basis vectorsNumber of atomsBCCFCCHCPSCDiamond1001021041061081010 0  100  200  300  400  500  600Number of basis vectorsNumber of atomsRocksaltZincblendeWurtziteFluoriteRutileC-perovskiteSpinel(a) (c)(b)Second order Third order Third orderFIG. 1. Number of basis vectors required to completely represent force constants in different crystal structures. (a) Forsecond-order force constants in elemental crystal structures: body-centered cubic (bcc), face-centered cubic (fcc), hexagonalclose-packed (hcp), simple cubic (sc), and diamond-type structures. (b) For third-order force constants in elemental crystalstructures. (c) For third-order force constants in ionic crystal structures, including rocksalt, zincblende, wurtzite, fluorite,rutile, cubic perovskite (C-perovskite), and spinel structures. The black solid lines indicate the theoretical values of 9N2 forsecond-order force constants and 27N3 for third-order force constants.1001021041061081010 0  200  400  600  800  1000Number of basis vectorsNumber of atomsNo cutoff distanceThird orderFIG. 2. Number of basis vectors to represent third-orderforce constants in the diamond-type structure, consideringvarious cutoff distances rc. The lattice constant of the dia-mond structure is denoted by a.1. Computational detailsSupercell force constants were determined by apply-ing standard linear regression to a displacement-forcedataset that consisted of Ns supercell structures. Theywere created by introducing random atomic displace-ments of magnitude d into the original supercell struc-ture with a given supercell size. The forces acting onatoms in these structures were then computed using theDFT calculation or polynomial MLPs [25, 26] represent-ing short-range interatomic interactions with a polyno-mial function of polynomial invariants for arbitrary rota-tions. In applications using DFT calculations, force cal-culations were performed using the plane-wave-basis pro-jector augmented wave method [27] within the Perdew–Burke–Ernzerhof exchange-correlation functional [28] asimplemented in the vasp code [29–31]. The other com-putational conditions are given in each application.The LTCs were calculated as the solution of thePeierls–Boltzmann equation under the relaxation timeapproximation [32–34] using the phono3py code [35, 36].To simplify the LTC calculations and measure the pre-cision of force constants, only phonon-phonon scatteringwas considered to obtain the phonon relaxation times.The phonon lifetimes were calculated from supercell forceconstants as the reciprocals of the imaginary part of thephonon self-energy corresponding to the bubble diagram,which were used as the relaxation times. Phonon prop-erties required for the LTC calculations such as phonongroup velocities and mode heat capacities were obtainedfrom dynamical matrices derived from second-order su-percell force constants. Mesh grids for sampling the re-ciprocal spaces are shown in each application.2. Estimation from DFT datasetsFirstly, we demonstrate an application of our ap-proach to diamond silicon, where displacement-forcedatasets were generated exclusively using DFT calcula-tions. In this application, we computed the forces actingon atoms in supercell structures with atomic displace-ments through DFT calculations and used them to esti-mate the force constants. The supercell structures wereconstructed using the 2 × 2 × 2 expansion of the unitcell. The cutoff energy was set to 300 eV, and the to-tal energies converged to less than 10−3 meV/supercell.12 0 25 50 75 100 125 150 0  10  20  30  40  50LTC (W/mK, 300 K)Number of supercells0.003 Å0.01 Å0.03 Å0.1 Å 0  200  400  600  800  1000Number of supercells(a) (b)FIG. 3. (a) Relationship between the LTC and the num-ber of supercells with displacement magnitudes d = 0.003,0.01, 0.03, and 0.1 Å in diamond silicon. Displacement-forcedatasets are derived exclusively from DFT calculations. Theminimum number of supercells needed to construct a full-rankpredictor matrix is five for the 2× 2× 2 supercell expansion.LTC values are shown for up to 50 supercells. (b) Relation-ship between LTC and the number of supercells with the samedisplacement magnitudes in diamond silicon, with LTC val-ues shown for up to 1000 supercells.We examined atomic displacements with magnitudes ofd = 0.003, 0.01, 0.03, and 0.1 Å.In Fig. 3, we present the relationship between theLTC, the magnitude of atomic displacements, and thenumber of supercells in the displacement-force dataset.When atomic displacements are larger than 0.01 Å andthe number of supercells is small, the LTC values are un-derestimated due to errors in the forces resulting from thetruncation of higher-order contributions. Consequently,minimal supercell structures are insufficient for accurateLTC prediction. On the other hand, accurate LTC valuescan be obtained even using a small dataset with atomicdisplacements of 0.003 Å. In the diamond structure, theatomic coordinates have no degrees of freedom, whichmeans that no residual forces arise. Consequently, nu-merical and truncation errors in the dataset are minimal.Then, we demonstrate LTC calculations using the DFTcalculation for wurzite-type ZnS. The wurtzite struc-ture affords a degrees of freedom in atomic coordinates,thus suggesting that the equilibrium structure optimizedthrough DFT calculations may possess residual forces ap-proaching the specified convergence tolerance. To createa displacement-force dataset, the forces acting on atomsin supercell structures were computed using only DFTcalculations. Supercell structures were constructed usingthe 3×3×2 expansion of the unit cell. The cutoff energywas set to 500 eV, and the total energies converged to lessthan 10−3 meV/supercell. The atomic positions in theunit cell were optimized until the residual forces were lessthan 0.0015 eV/Å. To examine the dependency of LTCon atomic displacements, we investigated magnitudes ofd = 0.003, 0.01, 0.03, and 0.1 Å.Figure 4 illustrates the relationship among LTC, thenumber of supercells, and the magnitude of atomic dis-placements in wurtzite-type ZnS. LTC values derived 0 25 50 75LTC (W/mK, 300 K)d = 0.003 Å d = 0.01 Å 0 25 50 75 0  250  500  750  1000LTC (W/mK, 300 K)Number of supercellsd = 0.03 Å 0  250  500  750  1000Number of supercellsd = 0.1 ÅResidual forces (eliminated) Residual forces (not eliminated) FIG. 4. Convergence behavior of LTC with respect to thenumber of supercells in the wurtzite-type ZnS. Displacement-force datasets are derived solely from DFT calculations. Thesupercell size is 3 × 3 × 2, with a minimum of 36 supercellsrequired to construct a full-rank predictor matrix. Averagevalues of the diagonal elements in LTC tensors are presented.The closed circles represent LTC values calculated using forceswithout subtracting residual forces, while the open circles rep-resent LTC values calculated using forces with residual forcesremoved.from raw displacement-force datasets are compared withthose derived from datasets where residual forces are sub-tracted. Residual forces can introduce numerical errors inforce constant calculations, significantly impacting LTCpredictions, especially when atomic displacements andresulting forces are small. For atomic displacements ofd = 0.003 and 0.01 Å, LTC values exhibit sensitivity towhether residual forces are subtracted from DFT forcesin supercell structures, as depicted in Fig. 4. Conversely,for displacements of d = 0.03 and 0.1 Å, the forces in su-percell structures are sufficiently large that LTC valuesappear independent of residual force subtraction.The results depicted in Fig. 4 highlight the importanceof subtracting residual forces for accurately and robustlyestimating force constants and LTCs. However, a notice-able discrepancy is observed in the converged LTC valueswhen using d = 0.003 Å compared to d = 0.01 and 0.03Å. This difference likely stems from residual forces andother numerical errors in the DFT calculations being rel-atively large compared to the magnitude of the forces.Conversely, for d = 0.1 Å, the significant contributionsfrom higher-order effects lead to slow LTC convergencewith respect to the number of supercells due to theirtruncation. Based on these findings, it is recommendedto use atomic displacements in the range of d = 0.01 to0.03 Å to obtain reliable LTC values in this case.In comparison to the diamond structure, the wurtzite13structure requires a greater number of supercell struc-tures due to the increased number of basis vectors forthe force constants. In such cases, adopting an efficientapproach for determining force constants proves highlyadvantageous [37]. This approach combines the cur-rent projector-based method with on-the-fly polynomialMLPs developed from DFT displacement-force datasets.DFT calculations are necessary only for supercell struc-tures used in MLP development, while force calculationsfor estimating force constants utilize the MLPs. More-over, even with small displacements in force calcula-tions, accurate force constants and LTC values can beachieved owing to the smooth potential energy surface ofMLPs. Consequently, this approach significantly reducesthe number of supercell structures needed for DFT calcu-lations, while maintaining accuracy comparable to forceconstant and LTC calculations obtained from DFT cal-culations for a much larger set of structures.3. Estimation from MLP datasetsHere, we perform force constant and LTC calcula-tions using displacement-force datasets obtained fromMLPs. The computational costs associated with MLPforce calculations are negligible in force constant andLTC calculations. We examine the performance of theprojector-based procedure for applications that are im-practical using only DFT calculations. MLPs offer sev-eral advantages: (1) They enable the generation of ex-tensive displacement-force datasets with significantly re-duced computational demands. (2) They exihibit smoothpotential energy surfaces around equilibrium structures,which allows us to employ small magnitudes of atomicdisplacements deriving supercell structures. (3) They fa-cilitate the elimination of residual forces in equilibriumstructures through local structure optimization.Force constant and LTC calculations were firstly con-ducted using MLP displacement-force datasets in dia-mond silicon. The polynomial MLP utilized in this studywas adopted from prior research [38], which investigatedLTCs at silicon grain boundaries through molecular dy-namics simulations. This MLP is also accessible in thePolynomial Machine Learning Potential Repository [39].It was trained on a dataset comprising diverse crystalstructures and their derivatives, incorporating randomatomic displacements and changes in cell shapes. ForLTC calculations, reciprocal spaces were sampled usinga 19× 19× 19 mesh grid.The current projector-based approach enables efficientestimation of force constants and LTCs for supercells con-taining up to 512 atoms, when no cutoff distance is ap-plied. Figure 5 illustrates the relationship between theLTC, supercell size, displacement magnitude, and thenumber of supercell structures within the displacement-force dataset. The LTC values obtained from the polyno-mial MLP using the 2×2×2 supercell are slightly smallerthan those obtained from the DFT calculation, as alsodemonstrated in Ref. 38. The complete third-order basisset sizes for the 2×2×2, 3×3×3, and 4×4×4 supercellexpansions are 777, 8800, and 49301, respectively. Con-sequently, the minimum numbers of supercells requiredto construct a full-rank predictor matrix X are 5, 14,and 33, respectively. When considering small displace-ment values, it is apparent that even a limited numberof supercells can yield accurate LTCs. Specifically, 5 and15 supercell structures generated from random displace-ments with d = 0.001 Å are adequate to obtain preciseLTCs for the 2× 2× 2 and 3× 3× 3 expansions, respec-tively. In contrast, minimal structures fail to provideaccurate LTC predictions when larger displacements areused. Due to the truncation of significant higher-ordercontributions, the LTC converges slowly as the numberof supercells increases, as evidenced in Fig. 5 (b).Figure 6 shows the relationship between the LTC andthe magnitude of atomic displacements. The LTC valuesare computed from displacement-force datasets compris-ing 10, 100, 1000, and 10000 supercell structures. Forvery small displacement values, such as d = 10−6 Å, theresulting forces are exceedingly small, which means thatnumerical noise can significantly affect the accuracy offorce constants. Consequently, a large number of super-cells are required to achieve accurate LTCs. To accu-rately evaluate third-order force constants without theneed for renormalizing higher-order contributions, it ispreferable to use atomic displacements within the rangeof 10−5 to 10−3 Å. Conversely, when atomic displace-ments exceed 0.01 Å, a substantial number of supercellsare also needed to achieve converged LTC values. Thisrequirement arises from truncation errors associated withhigher-order contributions, which complicate the compu-tation of converged LTCs.It is important to note that in this study, we utilized aconstant magnitude of atomic displacements to generateeach displacement-force dataset. However, the currentdiscussion is pertinent to displacement-force datasetsthat include a range of different magnitudes of atomicdisplacements. Such datasets are typically encounteredwhen estimating temperature-dependent force constantsrelevant to self-consistent phonons (e.g. [9, 10]), whichinherently include higher-order contributions to the forceconstants. Consequently, when evaluating self-consistentphonons, particularly at high temperatures, a large num-ber of supercell structures will be required.As described above, the imprecise prediction of LTCsis caused by errors in the forces and is closely relatedto the linear dependence among predictor variables inmatrix X. This phenomenon, known as multicollinear-ity, can render regression coefficients highly sensitive toerrors in the forces. In standard linear regression, thecondition number of matrix X⊤X is commonly used asan indicator of multicollinearity. The condition numberis determined by dividing the maximum eigenvalue ofX⊤X by the minimum eigenvalue. A significantly largecondition number suggests the presence of multicollinear-ity in the regression analysis.14 0 25 50 75 100 125 0  20  40  60  80  100LTC (W/mK, 300 K)d = 0.001 Å2x2x23x3x34x4x4 0 25 50 75 100 125 0  20  40  60  80  100Number of supercells (MLP)d = 0.003 Å2x2x23x3x34x4x4 0 25 50 75 100 125 0  20  40  60  80  100d = 0.01 Å2x2x23x3x34x4x4 0 25 50 75 100 125 0  1000  2000LTC (W/mK, 300 K)d = 0.03 Å2x2x23x3x34x4x4 0 25 50 75 100 125 0  1000  2000Number of supercells (MLP)d = 0.05 Å2x2x23x3x34x4x4 0 25 50 75 100 125 0  1000  2000d = 0.1 Å2x2x23x3x34x4x4(a)(b)FIG. 5. (a) Relationship between the LTC and the number of supercells with displacement magnitudes d = 0.001, 0.003,and 0.01 Å in diamond silicon. The minimum numbers of supercells required to construct a full-rank predictor matrix are 5,14, and 33 for the 2 × 2 × 2, 3 × 3 × 3, and 4 × 4 × 4 supercell expansions, respectively. LTC values are shown for up to 100supercells. (b) Relationship between LTC and the number of supercells with displacement magnitudes d = 0.03, 0.05, and 0.1Å in diamond silicon. LTC values are shown for up to 2000 supercells. 0 25 50 75 100 12510-6 10-5 10-4 10-3 10-2 10-1LTC (W/mK, 300 K)Displacement (Å)FIG. 6. Relationship between the LTC value and the mag-nitude of atomic displacements in diamond silicon with a su-percell size of 2 × 2 × 2. LTC values are calculated usingdisplacement-force datasets consisting of 10, 100, 1000, and10000 supercell structures.Figure 7 illustrates the dependence of the conditionnumber of matrix X⊤X for the third-order force con-stants with the 2 × 2 × 2 supercell on the number ofsupercell structures. The condition number is high whenusing only a small number of supercells, close to the min-imum number required to construct a full-rank predictormatrix. In such cases, the LTC predictions are highlysensitive to force errors, leading to imprecise LTC esti-mates. Conversely, the condition number decreases as thenumber of supercells increases, and this decrease alignswith the convergence patterns observed in LTCs. There-fore, using an orthonormal basis set for force constantsalong with a sufficiently large number of random displace-ments for force calculations creates optimal conditionsfor estimating force constants through standard linearregression.Next, we calculate LTCs in wurtzite AgI using a poly-nomial MLP, which was developed from a displacement-force DFT dataset using the on-the-fly approach [37].This dataset was made up of 300 supercell structureswith atomic displacements of 0.03 Å. We assess the forcesacting on atoms for Ns structures using the polynomialMLP. Here, the atomic positions optimized by the DFTcalculation are used as the equilibrium atomic positions.Therefore, residual forces at the equilibrium atomic posi-tions are eliminated from the forces to stabilize LTC cal-culations, particularly with small atomic displacements.15 0 10 20 30 40 0  200  400  600  800  1000Condition numberNumber of supercellsFIG. 7. Dependence of the condition number of X⊤X forthird-order force constants on the number of supercells in dia-mond silicon. The supercell size of 2×2×2 and the magnitudeof atomic displacements of d = 0.001 Å are used. The condi-tion number is nearly independent of the magnitude of atomicdisplacements, with similar values observed for d = 0.1 Å. 0 1 2 3 0  1000  2000LTC (W/mK, 300 K)Number of supercells (MLP)0.03 Å0.01 Å0.003 Å0.001 ÅFIG. 8. Convergence behavior of LTC with respect to thenumber of supercells in wurtzite AgI. The supercell size is 3×3×2, and the number of basis vectors for the third-order forceconstants is 7752, necessitating a minimum of 36 supercellstructures. The average values of the diagonal elements inLTC tensors are presented.Additionally, we sample the reciprocal spaces for LTCcalculations using the 19× 19× 10 mesh.Figure 8 shows the relationship between the LTC andthe number of supercells in wurtzite AgI. The supercellsare constructed using the 3× 3× 2 expansion of the unitcell. To construct a full-rank predictor matrix X, a min-imum of 36 supercell structures is required because thereare 7752 basis vectors for the third-order force constants.The convergence behavior of the LTC is influenced by themagnitude of atomic displacements. For example, a smallnumber of supercell structures are sufficient to achievea converged LTC value for small atomic displacements,analogous to the situation observed with diamond silicon.Conversely, for larger displacements, a greater number of 0 40 80 120 160 4  6  8  10  12LTC (W/mK, 300 K)Cutoff distance (Å)2x2x23x3x34x4x4FIG. 9. Dependence of the LTC on the cutoff distancein diamond silicon. The cutoff distances are applied only tothe third-order force constants. Force calculations for super-cell structures are performed using a polynomial MLP witha cutoff radius of 6 Å. The displacement magnitude is set tod = 0.001 Å and 100 supercell structures are used.supercell structures are necessary to obtain a precise LTCvalue that effectively incorporates higher-order contribu-tions.4. Estimation from MLP datasets using cutoff distanceFinally, we investigate the dependence of the LTC onthe cutoff distance for third-order force constants in dia-mond silicon. No cutoff distance is applied to the second-order force constants. Figure 9 shows how the LTC indiamond silicon varies with the cutoff distance. Forcecalculations for supercell structures are performed usingthe polynomial MLP described in the previous section.The displacement magnitude is set to d = 0.001 Å andthe number of supercell structures is 100. Although thepolynomial MLP has a cutoff radius of 6 Å, many-bodyterms of the polynomial MLP contribute to LTC predic-tions with cutoff distances exceeding 6 Å. For diamondsilicon, a cutoff distance of 8 Å in force constant calcula-tions yields reasonable LTC values. Thus, the introduc-tion of the cutoff distance significantly reduces the sizeof the basis set and computational demands while stillproviding accurate force constant calculations.VI. CONCLUSIONIn this study, we have proposed a projector-based ap-proach for efficiently constructing a complete orthonor-mal basis set for force constants. We have introduced sev-eral compression techniques to reduce the basis set sizeefficiently. These techniques involve compressing a pro-16jection matrix using eigenvectors of another projectionmatrix onto subspaces defined by partial requirementsfor the force constants. Our efficient procedures for basisset construction and force constant estimation, up to thefourth order, are implemented in the symfc code [16].The obtained basis set meets necessary requirementsfor the force constants, including permutation symme-try rules, sum rules, and symmetric properties dictatedby the space group operations of the target structure.Consequently, this complete orthonormal basis set facili-tates the precise determination of force constants, ensur-ing that they exactly satisfy the required constraints. Asdemonstrated in the applications of this study, proper-ties related to the force constants, such as LTCs, can becalculated both accurately and efficiently using this basisset.ACKNOWLEDGMENTSThis work was supported by JSPS KAKENHI GrantNumbers JP22H01756, and JP19H05787, JP24K08021,and JP21K04632.Appendix A: Derivation of projection matrix ontovector spaceIn this section, we derive the orthogonal projection ma-trix onto Ker(C⊤) [17, 18]. We begin by deriving theprojection matrix onto the vector space spanned by C,denoted as PC . Assuming that the column vectors ofC are linearly independent, the projection of a vector xonto the vector space spanned by C is expressed asPCx = Cα, (A1)where the projected vector is represented as a linear com-bination of the column vectors of C, with the coefficientsdenoted by the vectorα. Since the residual vector x−Cαis orthogonal to all column vectors of C, we have the re-lationship:C⊤[x−Cα] = 0. (A2)This can be transformed to yield C⊤Cα = C⊤x, leadingto the expression α = (C⊤C)−1C⊤x. By substitutingthis expression into Eq. (A1), we derive the projectionmatrix onto the vector space spanned by C asPC = C(C⊤C)−1C⊤. (A3)Since Ker(C⊤) constitutes the complementary subspaceof the vector space spanned by C, the orthogonal projec-tion matrix onto Ker(C⊤) can be expressed asPKer(C⊤) = I −C(C⊤C)−1C⊤, (A4)where I denotes the identity matrix. When C is an or-thogonal matrix, the orthogonal projection matrix ontoKer(C⊤) can be represented asPKer(C⊤) = I −CC⊤. (A5)Appendix B: Proofs of commutation relationsbetween projection matricesIn this section, we provide a proof of the commutationrelations between projection matrices, specifically show-ing that [Pspg,Pv] = 0 and [Pspg,Pw] = 0. Although weonly present the proof for second-order force constants,the argument can be extended to third-order force con-stants as well. Consider two orthogonal projection ma-trices, P1 and P2. The transpose of the product P1P2is given by (P1P2)⊤ = P⊤2 P⊤1 = P2P1. Therefore, ifP1P2 is symmetric, i.e., P1P2 = (P1P2)⊤, then P1 andP2 commute. We will now demonstrate that the productof two projection matrices is symmetric.(1) Proof of [Pspg,Pv] = 0. Firstly, we will show thatPvPspg is symmetric, PvPspg = (PvPspg)⊤. Since theprojection matrix Pv is given by Pv = I − CvC⊤v asshown in Eq. (22), the elements of P cv = CvC⊤v areexpressed asP cv,(iαjβ,i′α′j′β′) =1Nδii′δαα′δββ′ , (B1)where δii′ denotes the Kronecker delta. Therefore, whendenoting the product of P cv and Pspg as P ′ = P cvPspg,its elements are written asP ′(iαjβ,i′α′j′β′) =∑kγlδP cv,(iαjβ,kγlδ)Pspg(kγlδ,i′α′j′β′)=1N∑lPspg(iαlβ,i′α′j′β′) (B2)andP ′(i′α′j′β′,iαjβ) =∑kγlδP cv,(i′α′j′β′,kγlδ)Pspg(kγlδ,iαjβ)=1N∑lPspg(i′α′lβ′,iαjβ)=1N∑lPspg(iαjβ,i′α′lβ′). (B3)The final equality is derived from the symmetric propertyof Pspg. Using the definition of Pspg, these elements areexpressed asP ′(iαjβ,i′α′j′β′)=1N |G|∑l∑ĝ∈GΓ(ĝ)(iα,i′α′)Γ(ĝ)(lβ,j′β′) (B4)andP ′(i′α′j′β′,iαjβ)=1N |G|∑l∑ĝ∈GΓ(ĝ)(iα,i′α′)Γ(ĝ)(jβ,lβ′), (B5)respectively. The summations over index l in Eqs. (B4)and (B5) must then be equal, i.e.,∑lΓ(ĝ)(lβ,j′β′) =∑lΓ(ĝ)(jβ,lβ′), (B6)17because Γ(ĝ)(lβ,j′β′) can be non-zero only for single atoml such that operation ĝ moves atom l to atom j′ andits value corresponds to the representation of ĝ for theCartesian components (β, β′). When we denote thesesummations as Γ′(ĝ)(β,β′), the elements of P ′ are writtenasP ′(iαjβ,i′α′j′β′) = P ′(i′α′j′β′,iαjβ)=1N |G|∑ĝ∈GΓ(ĝ)(iα,i′α′)Γ′(ĝ)(β,β′), (B7)which is equivalent to P cvPspg = (P cvPspg)⊤. Therefore,PvPspg = (PvPspg)⊤ can also be derived, and it has beenproved that [Pspg,Pv] = 0.(2) Proof of [Pspg,Pw] = 0. We demonstrate thatPwPspg is symmetric. The projection matrix Pw is rep-resented by Pw = BwB⊤w . Therefore, the elements ofmatrix P ′′ = PwPspg are written asP ′′(iαjβ,i′α′j′β′) =12[Pspg(iαjβ,i′α′j′β′) + Pspg(jβiα,i′α′j′β′)](B8)andP ′′(i′α′j′β′,iαjβ) =12[Pspg(i′α′j′β′,iαjβ) + Pspg(j′β′i′α′,iαjβ)].(B9)Then, the projection matrix Pspg is orthogonal and hasthe following symmetric properties ofPspg(iαjβ,i′α′j′β′) = Pspg(i′α′j′β′,iαjβ) (B10)andPspg(jβiα,i′α′j′β′) = Pspg(i′α′j′β′,jβiα)= Pspg(j′β′i′α′,iαjβ), (B11)where the final equality is derived from the fact that Pspgis expressed as the sum of the Kronecker products ofrepresentation Γ(R̂). After applying these equations, wecan conclude that P ′′ is symmetric asP ′′(iαjβ,i′α′j′β′) = P ′′(i′α′j′β′,iαjβ), (B12)which proves that [Pspg,Pw] = 0.Appendix C: Implementation of force-constantestimationIn this section, we present a procedure for efficientlycomputing X⊤X and X⊤y when solving the normalequations derived from a displacement-force dataset. Thebasis set of the full projection matrix is obtained from twodecomposed matrices, as specified in Eq. (71). In thiscontext, the matrices Bspg∩w and Bv↓(spg∩w) appearingin Eq. (71) are denoted by D and E, respectively. Con-sequently, the basis sets B(FC2) and B(FC3) are writtenbyB(FC2) = D(FC2)E(FC2)B(FC3) = D(FC3)E(FC3). (C1)Using these expressions, the force constants are repre-sented byΘiα,jβ =∑bc(FC2)b B(FC2)iα,jβ,b=∑bc(FC2)b[∑dD(FC2)iα,jβ,dE(FC2)d,b](C2)andΘiα,jβ,kγ =∑bc(FC3)b B(FC3)iα,jβ,kγ,b=∑bc(FC3)b[∑dD(FC3)iα,jβ,kγ,dE(FC3)d,b],(C3)where cb denotes the expansion coefficient for b-th basis.We introduce the variables z(s,FC2)iα,d and z(s,FC3)iα,d definedbyz(s,FC2)iα,d = −∑jβu(s)jβ D(FC2)iα,jβ,dz(s,FC3)iα,d = −12∑jβ∑kγu(s)jβ u(s)kγD(FC3)iα,jβ,kγ,d. (C4)The forces acting on atoms in the supercell structure scan then be expressed asf(s)iα =∑bc(FC2)b∑dz(s,FC2)iα,d E(FC2)d,b+∑bc(FC3)b∑dz(s,FC3)iα,d E(FC3)d,b . (C5)In matrix form, the forces are represented asf (s) = Z(s)Ec= X(s)c, (C6)whereZ(s) =[Z(s,FC2) Z(s,FC3)]=· · ·... · · ·... · · ·· · · z(s,FC2)iα,d · · · z(s,FC3)iα,d · · ·· · ·... · · ·... · · · , (C7)andE =[E(FC2) 00 E(FC3)]. (C8)For the displacement-force dataset, the matrix Z is ex-pressed asZ =Z(1)Z(2)Z(3)...Z(Ns) =Z(1,FC2) Z(1,FC3)Z(2,FC2) Z(2,FC3)Z(3,FC2) Z(3,FC3)......Z(Ns,FC2) Z(Ns,FC3)=[Z(FC2) Z(FC3)]. (C9)18Using these representations, X⊤X and X⊤y can berewritten asX⊤X = E⊤Z⊤ZE =[E(FC2)⊤Z(FC2)⊤Z(FC2)E(FC2) E(FC2)⊤Z(FC2)⊤Z(FC3)E(FC3)E(FC3)⊤Z(FC3)⊤Z(FC2)E(FC2) E(FC3)⊤Z(FC3)⊤Z(FC3)E(FC3)], (C10)andX⊤y = E⊤Z⊤y =[E(FC2)⊤Z(FC2)⊤yE(FC3)⊤Z(FC3)⊤y], (C11)respectively. The matrix products with forms of Z⊤Zand Z⊤y are computed efficiently by partitioning Z andy into subsets, and X⊤X and X⊤y are calculated usingZ⊤Z and Z⊤y.Appendix D: Computational efficiencyIn this section, we present benchmark results that il-lustrate the computational efficiency of generating force-constant basis sets and computing force constants fromdisplacement-force datasets. Figure 10 depicts theelapsed time required for computing the basis sets ofsecond-order and third-order supercell force constantsfor both elemental and ionic structures. The computa-tional cost primarily depends on the number of atoms inthe supercell, as the number of elements in the second-order and third-order force-constant matrices scales with9N2 and 27N3, respectively. Additionally, the computa-tional cost varies according to the symmetry of the tar-get structure. Structures with lower symmetry generallyneed higher computational costs compared to those withhigher symmetry. This is because lower-symmetry struc-tures lack the symmetric properties that could efficientlyreduce the size of the basis sets. Consequently, solvingthe compressed eigenvalue problem as described by Eq.(70) becomes more computationally demanding for thesestructures.We conducted a benchmark to estimate second-and third-order force constants from displacement-forcedatasets using a polynomial MLP in diamond silicon.Figure 11 (a) illustrates the elapsed time and peak mem-ory usage associated with estimating force constants forvarious subset sizes. For this benchmark, we utilized adataset consisting of 3000 structures with a supercell sizeof 3 × 3 × 3 in diamond silicon. As depicted in Fig. 11(a), the computational cost remains relatively constantfor subset sizes up to 400 structures, whereas the mem-ory requirements increase approximately proportionallywith the subset size. Consequently, we recommend thatthe subset size should not exceed 100 to optimize perfor-mance.Figure 11 (b) displays the elapsed time and peak mem-ory usage for estimating force constants with a fixed sub-set size of 100 but varying numbers of supercells. Thecomputational cost scales linearly with the number of su-percells in the dataset, while the memory usage remainsnearly constant, regardless of the number of supercells.This result demonstrates that partitioning X⊤X effec-tively mitigates memory consumption.In Fig. 11 (c), we present the elapsed time and peakmemory usage for estimating force constants with a fixedsubset size of 100 and a fixed number of supercells (3000).We used three different supercell sizes: 2×2×2, 3×3×3,and 4 × 4 × 4. The results indicate that both compu-tational cost and memory requirements are significantlyinfluenced by the number of atoms in the supercell.Figure 11 (d) illustrates the elapsed time and peakmemory usage as a function of the cutoff distance. Forthis analysis, we used a subset size of 100 and 3000 super-cells. Introducing a cutoff distance for the force constantshelps reduce computational demands.[1] M. T. Dove, Introduction to Lattice Dynamics (Cam-bridge University Press, 1993).[2] Duane C. Wallace, Thermodynamics of crystals (DoverPublications, 1998).[3] John Michael Ziman, Electrons and phonons: the the-ory of transport phenomena in solids (Oxford UniversityPress, 1960).[4] K. Parlinski, Z. Q. Li, and Y. Kawazoe, “First-principlesdetermination of the soft mode in cubic zro2,” Phys. Rev.Lett. 78, 4063–4066 (1997).[5] O. Hellman, P. Steneteg, I. A. Abrikosov, andS. I. Simak, “Temperature dependent effective potentialmethod for accurate free energy calculations of solids,”Phys. Rev. B 87, 104111 (2013).[6] Terumasa Tadano and Shinji Tsuneyuki, “First-principles lattice dynamics method for strongly anhar-monic crystals,” J. Phys. Soc. Jpn. 87, 041015 (2018).[7] F. Eriksson, E. Fransson, and P. Erhart, “The hiphivepackage for the extraction of high-order force constantsby machine learning,” Adv. Theory Simul. 2, 1800184(2019).[8] Atsushi Togo, “First-principles phonon calculations withhttp://dx.doi.org/ 10.1093/acprof:oso/9780198507796.001.0001http://dx.doi.org/ 10.1093/acprof:oso/9780198507796.001.0001http://dx.doi.org/10.1103/PhysRevLett.78.4063http://dx.doi.org/10.1103/PhysRevLett.78.4063http://dx.doi.org/10.1103/PhysRevB.87.104111http://dx.doi.org/10.7566/JPSJ.87.041015http://dx.doi.org/https://doi.org/10.1002/adts.201800184http://dx.doi.org/https://doi.org/10.1002/adts.2018001841910-1100101102103104105106 0  300  600  900Elapsed time (s)Number of atomsBCCFCCHCPSCDiamond10-1100101102103104105106 0  200  400  600Elapsed time (s)Number of atomsRocksaltZincblendeWurtziteFluoriteRutileC-perovskiteSpinelFIG. 10. Elapsed time required for computing basis sets of second-order and third-order supercell force constants for elementaland ionic structures. Elapsed times were measured using the Intel(R) Xeon(R) Gold 6240 CPU. 0 600 1200 1800 2400 0  100  200  300  400 0 16 32 48 64Elapsed time (s)Peak memory (GB)Batch size 0 600 1200 1800 2400 0  2500  5000  7500 10000 0 16 32 48 64Elapsed time (s)Peak memory (GB)Number of supercells101102103104105106 0  200  400  600 0 64 128 192 256Elapsed time (s)Peak memory (GB)Number of atomsElapsed time (s)Peak memory (GB)Cutoff distance (Å)(a) (b)(d)(c)Peak memoryPeak memoryPeak memoryTimeTimeTime0600120018002400 4  6  8  10  12 0 16 32 48 64Peak memoryTimeFIG. 11. (a) Elapsed time and peak memory usage for estimating second- and third-order force constants with varyingsubset sizes in diamond silicon. The number of supercells is 3000, with a supercell size of 3 × 3 × 3 (216 atoms), and atomicdisplacements are set to d = 0.001 Å. (b) Elapsed time and peak memory usage for estimating second- and third-order forceconstants with a fixed subset size of 100 and varying numbers of supercells in diamond silicon. The supercell size is 3×3×3 (216atoms), and the magnitude of atomic displacements is d = 0.001 Å. (c) Elapsed time and peak memory usage for estimatingsecond- and third-order force constants with a fixed subset size of 100 and a fixed number of supercells (3000) in diamondsilicon. The magnitude of atomic displacements is d = 0.001 Å. (d) Dependence of elapsed time and peak memory usage onthe cutoff distance for estimating second- and third-order force constants with a fixed subset size of 100, a fixed number ofsupercells (3000), and atomic displacements of d = 0.001 Å. Elapsed times were measured using the Intel(R) Xeon(R) Gold6240 CPU.20phonopy and phono3py,” J. Phys. Soc. Jpn. 92, 012001(2023).[9] Ion Errea, Matteo Calandra, and Francesco Mauri,“First-principles theory of anharmonicity and the inverseisotope effect in superconducting palladium-hydride com-pounds,” Phys. Rev. Lett. 111, 177002 (2013).[10] Ambroise van Roekeghem, Jesús Carrete, and NatalioMingo, “Quantum self-consistent ab-initio lattice dynam-ics,” Comput. Phys. Commun. 263, 107945 (2021).[11] J. F. Nye, Physical Properties of Crystals (Oxford Uni-versity Press, 1985).[12] M. El-Batanouny and F. Wooten, Symmetry and Con-densed Matter Physics: A Computational Approach(Cambridge University Press, 2008).[13] G. Venkataraman, L. A. Feldkamp, and V. C. Sahni,Dynamics of perfect crystals (MIT press, 1975).[14] Fei Zhou, Weston Nielson, Yi Xia, and VidvudsOzoliņš, “Lattice anharmonicity and thermal conductiv-ity from compressive sensing of first-principles calcula-tions,” Phys. Rev. Lett. 113, 185501 (2014).[15] T. Tadano and S. Tsuneyuki, “Self-consistent phonon cal-culations of lattice dynamical properties in cubic srtio3with first-principles anharmonic force constants,” Phys.Rev. B 92, 054301 (2015).[16] Atsuto Seko and Atsushi Togo, “Symfc,” https://github.com/symfc/symfc.[17] Gilbert Strang, Introduction to Linear Algebra, 6th ed.(CUP, 2023).[18] Haruo Yanai, Kei Takeuchi, and Yoshio Takane, Projec-tion Matrices, Generalized Inverse Matrices, and Singu-lar Value Decomposition, Statistics for social and behav-ioral sciences (Springer, Dordrecht, 2011).[19] R. Piziak, P.L. Odell, and R. Hahn, “Constructingprojections on sums and intersections,” Comput. Math.Appl. 37, 67–74 (1999).[20] John von Neumann, “On rings of operators. reductiontheory,” Annals of Math. 50, 401 (1949).[21] Israel Halperin, “The product of projection operators,”Acta Sci. Math.(Szeged) 23, 96–99 (1962).[22] Aurél Galántai, Projectors and projection methods, Vol. 6(Springer Science & Business Media, 2013).[23] Robert Sedgewick, Algorithms in c, part 5: graph al-gorithms, third edition, 3rd ed. (Addison-Wesley Profes-sional, 2001).[24] G. Strang, Linear Algebra and Learning from Data(Wellesley-Cambridge Press, 2019).[25] Atsuto Seko, Atsushi Togo, and Isao Tanaka, “Group-theoretical high-order rotational invariants for structuralrepresentations: Application to linearized machine learn-ing interatomic potential,” Phys. Rev. B 99, 214108(2019).[26] Atsuto Seko, “Tutorial: Systematic development of poly-nomial machine learning potentials for elemental and al-loy systems,” J. Appl. Phys. 133, 011101 (2023).[27] P E Blöchl, “Projector augmented-wave method,” Phys.Rev. B 50, 17953–17979 (1994).[28] J P Perdew, K Burke, and M Ernzerhof, “Generalizedgradient approximation made simple,” Phys. Rev. Lett.77, 3865–3868 (1996).[29] G Kresse and J Hafner, “Ab initio molecular-dynamicsfor liquid-metals,” Phys. Rev. B 47, 558–561 (1993).[30] G Kresse and J Furthmüller, “Efficient iterative schemesfor ab initio total-energy calculations using a plane-wavebasis set,” Phys. Rev. B 54, 11169–11186 (1996).[31] G Kresse and D Joubert, “From ultrasoft pseudopoten-tials to the projector augmented-wave method,” Phys.Rev. B 59, 1758–1775 (1999).[32] R. E. Peierls, “Zur kinetischen theorie der wärmeleitungin kristallen,” Ann. Phys. 395, 1055–1101 (1929).[33] R. E. Peierls, Quantum theory of solids (Oxford Univer-sity Press, 2001).[34] Philip B. Allen and Vasili Perebeinos, “Temperature ina peierls-boltzmann treatment of nonlocal phonon heattransport,” Phys. Rev. B 98, 085427 (2018).[35] A. Togo, L. Chaput, and I. Tanaka, “Distributions ofphonon lifetimes in brillouin zones,” Phys. Rev. B 91,094306 (2015).[36] Atsushi Togo, Laurent Chaput, Terumasa Tadano, andIsao Tanaka, “Implementation strategies in phonopy andphono3py,” J. Phys. Condens. Matter 35, 353001 (2023).[37] Atsushi Togo and Atsuto Seko, “On-the-fly training ofpolynomial machine learning potentials in computing lat-tice thermal conductivity,” J. Chem. Phys. 160, 211001(2024).[38] Susumu Fujii and Atsuto Seko, “Structure and latticethermal conductivity of grain boundaries in silicon byusing machine learning potential and molecular dynam-ics,” Comput. Mater. Sci. 204, 111137 (2022).[39] A. Seko, Polynomial Machine Learning Potential Repos-itory at Kyoto University, https://sekocha.github.io.http://dx.doi.org/ 10.7566/JPSJ.92.012001http://dx.doi.org/ 10.7566/JPSJ.92.012001http://dx.doi.org/ 10.1103/PhysRevLett.111.177002http://dx.doi.org/ https://doi.org/10.1016/j.cpc.2021.107945http://dx.doi.org/10.1017/CBO9780511755736http://dx.doi.org/10.1017/CBO9780511755736http://dx.doi.org/10.1103/PhysRevLett.113.185501http://dx.doi.org/10.1103/PhysRevB.92.054301http://dx.doi.org/10.1103/PhysRevB.92.054301https://github.com/symfc/symfchttps://github.com/symfc/symfchttp://dx.doi.org/10.1007/978-1-4419-9887-3http://dx.doi.org/10.1007/978-1-4419-9887-3http://dx.doi.org/10.1007/978-1-4419-9887-3http://dx.doi.org/ https://doi.org/10.1016/S0898-1221(98)00242-9http://dx.doi.org/ https://doi.org/10.1016/S0898-1221(98)00242-9https://api.semanticscholar.org/CorpusID:124439084http://dx.doi.org/ 10.1103/PhysRevB.99.214108http://dx.doi.org/ 10.1103/PhysRevB.99.214108http://dx.doi.org/10.1063/5.0129045http://dx.doi.org/ https://doi.org/10.1002/andp.19293950803http://dx.doi.org/10.1103/PhysRevB.98.085427http://dx.doi.org/10.1088/1361-648X/acd831http://dx.doi.org/10.1063/5.0211296http://dx.doi.org/10.1063/5.0211296http://dx.doi.org/ https://doi.org/10.1016/j.commatsci.2021.111137https://sekocha.github.io Projector-based efficient estimation of force constants Abstract Introduction Projection matrices for force constants Crystal potential and force constants Projection matrix for second-order force constants Projection matrix for third-order force constants Efficient basis set construction Decomposition of projection matrix Non-commutative projection matrices Compression of projection matrices Application of sum rules Application of permutation symmetry rules Effective use of lattice translation group Properties of projection matrices Eigenvalue problems for projection matrices Procedure to calculate eigenvectors of full projection matrix Practical approximations Force constant estimation Force-constant basis set expansion Force constant estimation Results and discussion Basis set construction Complete basis set Application of cutoff distance Force constant estimation Computational details Estimation from DFT datasets Estimation from MLP datasets Estimation from MLP datasets using cutoff distance Conclusion Acknowledgments Derivation of projection matrix onto vector space Proofs of commutation relations between projection matrices Implementation of force-constant estimation Computational efficiency References