The small displacement method

We want to devise now a practical method to compute the phonon frequencies in a crystal, which we can then use to compute its free energy. In the previous section we have seen that these frequencies are the square roots of the eigenvalues of the dynamical matrix, which is the lattice Fourier transform of the force constant matrix. Therefore, knowledge of the force constant matrix is sufficient to compute the whole phonon spectrum. This can be easily computed if the potential energy function is know analytically. If it is not, but we are able to obtain forces by some other means, for example using computer simulation, then we can exploit the relation between ${\bf\Phi}$ and the forces, given by . Let us consider the crystal in the ground state and displace only one particle from its equilibrium position, say the one in the origin. Let the displacement be along the $x$-axis, for simplicity, ${\bf u} = (u,0, 0)$. The forces acting on all other particles in the system are given by , and they are:

$$(f_i^x,f_i^y,f_i^z) = - {\bf\Phi}({\bf r}^0_i) \cdot (u,0,0) = -u \begin{pmatrix} \phi_{0i}^{xx} \\ \phi_{0i}^{yx}\\ \phi_{0i}^{zx} \end{pmatrix}$$ (7.62)

from which we obtain:

$$\phi_{0i}^{xx} = - \frac{f_i^x}{u}; \quad \quad \phi_{0i}^{yx} = - \frac{f_i^y}{u}; \quad \quad \phi_{0i}^{zx} = - \frac{f_i^z}{u}.$$ (7.63)

This is the first column of the force constant matrix. To obtain the other two columns we simply need to repeat the procedure, but this time with displacements $(0,u,0)$ and $(0,0,u)$. A set of three forces calculations, one for each displacement, is therefore all that is needed to compute the full force constant matrix. However, in high symmetry systems the number of displacements can be reduced. For example, if the system has a symmetry that makes the $x$ and the $y$ axis equivalent, the forces induced by a displacement along the $y$-axis can simply be obtained by applying the symmetry operation to the forces induced by the first displacement along the $x$-axis. For cubic systems for example, only a single displacement is ever needed.

In a practical calculation it is not possible to include all (infinite) particles and therefore it is not possible to compute the forces on all particles and obtain the full force constant matrix. The usual approach is similar to the one described for the linear chain in the previous section, i.e. we consider a supercell, which is a multiple $i \times j \times k$ of the primitive cell. The periodicity is now described by new lattice vectors:

$${\bf A}_1 = i {\bf a}_1; \quad \quad {\bf A}_2 = j {\bf a}_2; \quad \quad {\bf A}_3 = k {\bf a}_3,$$ (7.64)

and reciprocal lattice vectors:

$${\bf B}_1 = 2\pi \frac{ {\bf A}_2 \times {\bf A}_3}{V}; \quad {\b... ...es {\bf A}_1}{V}; \quad {\bf B}_3 = 2\pi \frac{ {\bf A}_1 \times {\bf A}_2}{V},$$ (7.65)

with $V = ({\bf A}_1 \times {\bf A}_2) \cdot {\bf A}_3$ the volume of the supercell. This means that if we displace one particle by an amount ${\bf u}$ in the origin, we also displace all its periodic images at positions ${\bf L} = n {\bf A}_1 + m {\bf A}_2 + l {\bf A}_3$ with $n,m,l$ any three integers. These displacements are all equal to ${\bf u}$ and therefore the forces that we are obtaining are:

$${\bf f}_i = - \sum_{\bf L} {\bf\Phi}({\bf r}^0_i + {\bf L}) \cdot {\bf u},$$ (7.66)

where the sum runs over all possible superlattice vectors ${\bf L}$. It is clear, therefore, that we do not have access to the single elements of ${\bf\Phi}({\bf r}^0_i)$, but only the supercell force constant matrix, defined as:

$${\bf\Phi}^{SC}({\bf r}^0_i)= \sum_{\bf L} {\bf\Phi}({\bf r}^0_i + {\bf L}).$$ (7.67)

This also implies that it is impossible to obtain the exact dynamical matrix, which requires knowledge of every element of the force constant matrix. Let us rewrite  as:

$${\bf D}({\bf q}) = \frac{1}{m} \sum_j \sum_{\bf L} e^{i {\bf q}\cdot ({\bf r}^0_j + {\bf L})} {\bf\Phi}({\bf r}^0_j+{\bf L}) =$$

$$\frac{1}{m} \sum_j e^{i {\bf q}\cdot {\bf r}^0_j} \sum_{\bf L} e^{i {\bf q}\cdot {\bf L}} {\bf\Phi}({\bf r}^0_j+{\bf L})$$ (7.68)

with the sum over $j$ running over the lattice sites ${\bf r}^0_j$ belonging to the supercell. Now consider the term $\sum_{\bf L} e^{i {\bf q}\cdot {\bf L}} {\bf\Phi}({\bf r}^0_j+{\bf L})$ and compare with Eq. . For every ${\bf q}$ such that $e^{i {\bf q}\cdot {\bf L}} = 1$ we have

$${\bf D}({\bf q}) = \frac{1}{m} \sum_j e^{i {\bf q}\cdot {\bf r}^0_j} {\bf\Phi}^{SC}({\bf r}^0_j).$$ (7.69)

Since ${\bf\Phi}^{SC}$ can be calculated, at these particular ${\bf q}$ vectors the dynamical matrix is exact. These are the ${\bf q}$ vectors that are linear combinations of the reciprocal vectors of the supercell, because if ${\bf q} = n^\prime {\bf B}_1 + m^\prime {\bf B}_2 + l^\prime {\bf B}_3$ then for any ${\bf L} = n {\bf A}_1 + m {\bf A}_2 + l {\bf A}_3$ we have ${\bf q} \cdot {\bf L} = (n n^\prime + m m^\prime + l l^\prime)2 \pi$, which gives $e^{i {\bf q}\cdot {\bf L}} = 1$.

To illustrate this property, let us consider a simple cubic system, with lattice and reciprocal vectors given by:

$${\bf a}_1 = (1,0,0); \quad \quad {\bf a}_2 = (0,1,0); \quad \quad {\bf a}_3 = (0,0,1)$$

$$\frac{{\bf b}_1}{2\pi} = (1,0,0); \quad \quad \frac{{\bf b}_2}{2\pi} = (0,1,0); \quad \quad \frac{{\bf b}_3}{2\pi} = (0,0,1).$$ (7.70)

Let us build a $2 \times 2 \times 2$ supercell:

$${\bf A}_1 = 2 {\bf a}_1; \quad \quad {\bf A}_2 = 2 {\bf a}_2; \quad \quad {\bf A}_3 = 2 {\bf a}_3$$

$${\bf B}_1 = \frac{{\bf b}_1}{2}; \quad \quad {\bf B}_2 = \frac{{\bf b}_2}{2}; \quad \quad {\bf B}_3 = \frac{{\bf b}_3}{2}$$ (7.71)

With these settings, using the procedure described above we obtain exact phonon frequencies at ${\bf q} = \frac{{\bf b}_1}{2}$, for example, which is a zone boundary phonon. We also have an exact phonon at ${\bf q} = -\frac{{\bf b}_1}{2}$, which is another zone boundary phonon and all other combinations along the other axis.

If we made a $3 \times 3 \times 3$ supercell, then we would still have one wave-vector (for any cartesian direction) between zero (the $\Gamma$ point) and zone boundary at which the phonons are exact: ${\bf q} = \frac{{\bf b}_1}{3}$, but note that there would be no zone-boundary phonons. This is equivalent to the statement in  where we had the $N/2$ term only for $N$ even. Another way of saying this is that an 'odd' supercell cannot accommodate stationary modes.

For a $4 \times 4 \times 4$ supercell the exact phonons are at ${\bf q} = \pm \frac{{\bf b}_1}{4}; \pm \frac{{\bf b}_1}{2}$ and so on. We recover the zone boundary phonons and, in addition, we get exact phonons also in between $\Gamma$ and zone boundary. Of course, every supercell, including the primitive cell, has exact (and trivial) phonons at $\Gamma$.

The property described above is quite useful, because it allows for an independent check on the consistency of the calculations. If we repeat them with both a $2 \times 2 \times 2$ and a $4 \times 4 \times 4$ supercells, the two sets of calculations must produce the same phonons at zone boundary.

We see that as we increase the size of the supercell we populate the BZ with more and more exact phonons and in between those points we obtain a Fourier interpolation, which becomes more and more accurate as we increase the size of the supercell. Eventually the supercell is so large that the force constant matrix is negligible at its edges. When this happens, the only term contributing in the sum in  is that with ${\bf L} = 0$, and ${\bf\Phi}({\bf r}^0_i) \approx {\bf\Phi}^{SC}({\bf r}^0_i)$, with the $\approx$ sign becoming an equal sign in the limit of the size of the supercell going to infinity. In this limit the Fourier interpolation is accurate everywhere. Note that the sum in  may include vectors on the edges of the Wigner-Seitz cell, and therefore only one of these terms has to be included in the sum. Alternatively, one can divide the value of ${\bf\Phi}^{SC}({\bf r}^0_i)$ by the number of equivalent vectors on the edges of the WS cell and include all of them in the sum .