Formalism for program alm

Interatomic force constants (IFCs)

The starting point of the computational methodology is to approximate the potential energy of interacting atoms by a Taylor expansion with respect to atomic displacements by

(1)\[\begin{split}&U - U_{0} = \sum_{n=1}^{N} U_{n} = U_{1} + U_{2} + U_{3} + \cdots, \\ &U_{n} = \frac{1}{n!} \sum_{\substack{\ell_{1}\kappa_{1}, \dots, \ell_{n}\kappa_{n} \\ \mu_{1},\dots, \mu_{n}}} \Phi_{\mu_{1}\dots\mu_{n}}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n}) \; u_{\mu_{1}}(\ell_{1}\kappa_{1})\cdots u_{\mu_{n}}(\ell_{n}\kappa_{n}).\end{split}\]

Here, \(u_{\mu}(\ell\kappa)\) is the atomic displacement of \(\kappa\)th atom in the \(\ell\)th unit cell along \(\mu\)th direction, and \(\Phi_{\mu_{1}\dots\mu_{n}}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n})\) is the \(n\)th-order interatomic force constant (IFC).

Symmetry relationship between IFCs

The are several relationships between IFCs which may be used to reduce the number of independence IFCs.

  • Permutation

Firstly, IFCs should be invariant under the exchange of triplet \((\ell,\kappa,\mu)\), e.g.

  • Periodicity

Secondly, since IFCs should depend on interatomic distances, they are invariant under a translation in units of lattice vector, namely

  • Crystal symmetry

A crystal symmetry operation maps an atom \(\vec{r}(\ell\kappa)\) to another equivalent atom \(\vec{r}(LK)\) by rotation and translation. Since the potential energy is invariant under any crystal symmetry operations, IFCs should transform under a symmetry operation as follows:

(2)\[\sum_{\nu_{1},\dots,\nu_{n}}\Phi_{\nu_{1}\dots\nu_{n}}(L_{1}K_{1};\dots;L_{n}K_{n}) O_{\mu_{1}\nu_{1}}\cdots O_{\mu_{n}\nu_{n}} = \Phi_{\mu_{1}\dots\mu_{n}}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n}),\]

where \(O\) is the rotational matrix of the symmetry operation. Let \(N_{s}\) be the number of symmetry operations, there are \(N_{s}\) relationships between IFCs which may be used to find independent IFCs.


In the current implementation of alm, independent IFCs are searched in Cartesian coordinate where the matrix element \(O_{\mu\nu}\) is 0 or \(\pm1\) in all symmetry operations except for those of hexagonal (trigonal) lattice. Also, except for hexagonal (trigonal) systems, the product \(O_{\mu_{1}\nu_{1}}\cdots O_{\mu_{n}\nu_{n}}\) in the left hand side of equation (2) becomes non-zero only for a specific pair of \(\{\nu\}\) (and becomes 0 for all other \(\{\nu\}\)s). Therefore, let \(\{\nu^{\prime}\}\) be such a pair of \(\{\nu\}\), the equation (2) can be reduced to

(3)\[\Phi_{\nu_{1}^{\prime}\dots\nu_{n}^{\prime}}(L_{1}K_{1};\dots;L_{n}K_{n}) = s \Phi_{\mu_{1}\dots\mu_{n}}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n}),\]

where \(s=\pm1\). The code employs equation (3) instead of equation (2) to reduce the number of IFCs. If IFCs of the left-hand side and the right-hand side of equation (3) are equivalent and the coupling coefficient is \(s=-1\), the IFC is removed since it becomes zero. For hexagonal (trigonal) systems, there can be symmetry operations where multiple terms in the left-hand side of equation (2) become non-zero. For such cases, equation (2) is not used to reduce the number of IFCs. Alternatively, the corresponding symmetry relationships are imposed as constraints between IFCs in solving fitting problems.

Constraints between IFCs

Since the potential energy is invariant under rigid translation and rotation, it may be necessary for IFCs to satisfy corresponding constraints.

The constraints for translational invariance are given by

(4)\[\sum_{\ell_{1}\kappa_{1}}\Phi_{\mu_{1}\mu_{2}\dots\mu_{n}}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2};\dots;\ell_{n}\kappa_{n}) = 0,\]

which should be satisfied for arbitrary pairs of \(\ell_{2}\kappa_{2},\dots,\ell_{n}\kappa_{n}\) and \(\mu_{1},\dots,\mu_{n}\). The code alm imposes equation (4) by default (ICONST = 1).

The constraints for rotational invariance are

\[\begin{split}&\sum_{\ell^{\prime}\kappa^{\prime}} (\Phi_{\mu_{1}\dots\mu_{n}\nu}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n};\ell^{\prime}\kappa^{\prime}) r_{\mu}(\ell^{\prime}\kappa^{\prime}) - \Phi_{\mu_{1}\dots\mu_{n}\mu}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n};\ell^{\prime}\kappa^{\prime}) r_{\nu}(\ell^{\prime}\kappa^{\prime})) \\ & \hspace{10mm} + \sum_{\lambda = 1}^{n}\sum_{\mu_{\lambda}^{\prime}} \Phi_{\mu_{1}\dots\mu_{\lambda}^{\prime}\dots\mu_{n}}(\ell_{1}\kappa_{1};\dots;\ell_{\lambda}\kappa_{\lambda};\dots;\ell_{n}\kappa_{n}) (\delta_{\mu,\mu_{\lambda}}\delta_{\nu,\mu_{\lambda}^{\prime}} - \delta_{\nu,\mu_{\lambda}}\delta_{\mu,\mu_{\lambda}^{\prime}}) = 0,\end{split}\]

which must be satisfied for arbitrary pairs of \((\ell_{1}\kappa_{1},\dots,\ell_{n}\kappa_{n};\mu_{1},\dots,\mu_{n};\mu,\nu)\). This is complicated since \((n+1)\)th-order IFCs (first line) are related to \(n\)th-order IFCs (second line).

For example, the constraints for rotational invariance related to harmonic terms can be found as

(5)\[\begin{split}&\sum_{\ell_{2}\kappa_{2}} (\Phi_{\mu_{1}\nu}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2})r_{\mu}(\ell_{2}\kappa_{2})-\Phi_{\mu_{1}\mu}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2})r_{\nu}(\ell_{2}\kappa_{2})) \notag \\ & \hspace{10mm} + \Phi_{\nu}(\ell_{1}\kappa_{1})\delta_{\mu,\mu_{1}} - \Phi_{\mu}(\ell_{1}\kappa_{1})\delta_{\nu,\mu_{1}} = 0,\end{split}\]


(6)\[\begin{split}& \sum_{\ell_{3}\kappa_{3}} (\Phi_{\mu_{1}\mu_{2}\nu}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2};\ell_{3}\kappa_{3}) r_{\mu}(\ell_{3}\kappa_{3}) \notag - \Phi_{\mu_{1}\mu_{2}\mu}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2};\ell_{3}\kappa_{3}) r_{\nu}(\ell_{3}\kappa_{3})) \\ & \hspace{10mm} + \Phi_{\nu\mu_{2}}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2})\delta_{\mu,\mu_{1}} - \Phi_{\mu\mu_{2}}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2})\delta_{\nu,\mu_{1}} \notag \\ & \hspace{10mm} + \Phi_{\mu_{1}\nu}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2})\delta_{\mu,\mu_{2}} - \Phi_{\mu_{1}\mu}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2})\delta_{\nu,\mu_{2}} = 0.\end{split}\]

When NORDER = 1, equation (5) will be considered if ICONST = 2, whereas equation (6) will be neglected. To further consider equation (6), please use ICONST = 3, though it may enforce a number of harmonic IFCs to be zero since cubic terms don’t exist in harmonic calculations (NORDER = 1).

Estimate IFCs by least-square fitting

The code alm extracts harmonic and anharmonic IFCs from a displacement-force data set by solving the following linear least-square problem:

(7)\[\text{minimize} \ \ \chi^{2} = \sum_{t}^{m} \sum_{i} \|F_{i,t}^{\mathrm{DFT}} - F_{i,t}^{\mathrm{ALM}} \|^{2}.\]

Here, \(m\) is the number of atomic configurations and the index \(i = (\ell,\kappa,\mu)\) is the triplet of coordinates. The model force \(F_{i,t}^{\mathrm{ALM}}\) is a linear function of IFCs \(\{\Phi\}\) which can be obtained by differentiating \(U\) (Eq. (1)) by \(u_{i}\). The parameters (IFCs) are determined so as to best mimic the atomic forces obtained by DFT calculations, \(F_{i,t}^{\mathrm{DFT}}\).

To evaluate goodness of fit, alm reports the relative error \(\sigma\) defined by

(8)\[\sigma = \sqrt{\frac{\chi^{2}}{\sum_{t}^{m}\sum_{i} (F_{i,t}^{\mathrm{DFT}})^{2}}},\]

where the numerator is the residual of fit and the denominator is the square sum of DFT forces.