Generation of force constants

A displacement is represented as $3 \times 9$ matrix for easy treatment as shown in below,

$$\mathbf{U} = \begin{pmatrix}U_{\alpha} & U_{\beta} & U_{\gamma} &... ...\ 0 & 0 & 0 & 0 & 0 & 0 & U_{\alpha} & U_{\beta} & U_{\gamma} \\ \end{pmatrix}.$$ (3.3)

First we must know the site point symmetry of the atom which is displaced. To symmetrize force constants, similarity transformation  [12] is used like,

$$\mathbf{P}^{\mathrm{transform}} = \mathbf{S} \cdot \mathbf{P}^{\mathrm{original}} \cdot \mathbf{S}^{-1},$$ (3.4)

where $\mathbf{S}$ is a site point symmetry operation where the transformation center is on the displaced atom and is $3 \times 3$ matrix. In this equation, force constants are treated as $3 \times 3$ matrix. After transformation, force constants are converted to $9 \times 1$ matrix. The position of atom $\mathbf{r}(jl)$ is also translated to $\mathbf{r}^{\mathrm{trans}} = \mathbf{S \cdot r}(jl)$. These matrices are connected in one matrix equation as follows,

$$\begin{multline} \mathbf{F} \begin{pmatrix} F_{\alpha}^{1}(\mathbf{r}^{1}) \\... ...ma\gamma}^{n_\mathrm{op}} \\ \end{pmatrix} \\ \end{pmatrix}, \end{multline}$$

where, the sufficient number of directions where the atom is displaced are denoted as $n_{\mathrm{dir}}$ and the number of the site point symmetry operations are denoted as $n_{\mathrm{op}}$. At this point, force constants are unknown, but symmetry operations are known. In this case, the symmetry operation should be done to displacement matrices rather than to force constants. Therefore the symmetry operations are separated from transformed force constants and the separated symmetry properties are connected to the displacement matrices as follows,

$$\mathbf{F} = - (\mathbf{U} \ast \mathbf{A}) \cdot \mathbf{P},$$ (3.5)

where $\mathbf{A}$ is the symmetry property matrix which is determined by rule. The operation '$\ast$' represents the operation between $\mathbf{U}$ and $\mathbf{A}$. This equation is described as,

$$\mathbf{F} \begin{pmatrix}\mathbf{F}^1(\mathbf{r}^1) \\ \vdots \\... ...a} \\ P_{\gamma\alpha} \\ P_{\gamma\beta} \\ P_{\gamma\gamma} \\ \end{pmatrix},$$ (3.6)

where $\mathbf{A}^{n}$ is $9 \times 9$ symmetry property matrix.

This equation can be solved by matrix inversion, as shown below,

$$\mathbf{P} = - (\mathbf{U} \ast \mathbf{A})^{-1} \cdot \mathbf{F}.$$ (3.7)

In most cases, force constants $\mathbf{P}$ are over-determined. In fropho, the least-squares solution of force constants is calculated by using singular value decomposition [13] in LAPACK routine [16].