Calculation of the force constant matrix

We calculate $\Phi_{i s \alpha , j t \beta}$ by the small-displacement method. In harmonic approximation the $\alpha$ Cartesian component of the force exerted on the atom at position ${\bf R}_i^0 + {\bf\tau}_s$is

$$F_{i s \alpha} = - \sum_{j t \beta} \Phi_{i s \alpha , j t \beta}~u_{j t \beta}$$ (2)

where $u_{j s \beta}$ is the displacement of the atom in ${\bf R}_j^0 + {\bf\tau}_t$ along the direction $\beta$. The force constant matrix can be calculated by

$$\Phi_{i s \alpha, j t \beta}=-\frac{F_{i s \alpha, j t \beta}}{ u_{j t \beta} }$$ (3)

by displacing once at a time all the atoms of the lattice along the three Cartesian components by $u_{j t \beta}$, and calculating the forces $F_{i s \alpha, j t \beta}$ induced on the atoms in ${\bf R}_i^0 + {\bf\tau}_s$. Since the crystal is invariant under traslations of any lattice vector, it is only necessary to displace the atoms in one primitive cell and calculate the forces induced on all the other atoms of the crystal. In what follows we will assume this as understood and put simply j = 0.

It is important to appreciate that the $\Phi_{l s \alpha , l^\prime t \beta}$ in the formula for $D_{s \alpha , t \beta} ( {\bf k} )$ is the force-constant matrix in the infinite lattice, with no restriction on the wavevector ${\bf k}$, whereas the calculations of $\Phi_{l s \alpha , l^\prime t \beta}$ can only be done in supercell geometry. Without a further assumption, it is strictly impossible to extract the infinite-lattice $\Phi_{l s \alpha , l^\prime t \beta}$from supercell calculations, since the latter deliver information only at wavevectors that are reciprocal lattice vectors of the superlattice (see note below ¹). The further assumption needed is that the infinite-lattice $\Phi_{l s \alpha , l^\prime t \beta}$ vanishes when the separation ${\bf R}_{l^\prime t} - {\bf R}_{l s}$ is such that the positions ${\bf R}_{l s}$ and ${\bf R}_{l^\prime t}$ lie in different Wigner-Seitz (WS) cells of the chosen superlattice. More precisely, if we take the WS cell centred on ${\bf R}_{l^\prime t}$, then the infinite-lattice value of $\Phi_{l s \alpha , l^\prime t \beta}$ vanishes if ${\bf R}_{l s}$ is in a different WS cell; it is equal to the supercell value if ${\bf R}_{l s}$ is wholly within the same WS cell; and it is equal to the supercell value divided by an integer P if ${\bf R}_{l s}$ lies on the boundary of the same WS cell, where P is the number of WS cells having ${\bf R}_{l s}$ on their boundary. With this assumption, the $\Phi_{l s \alpha , l^\prime t \beta}$ elements will converge to the correct infinite-lattice values as the dimensions of the supercell are systematically increased.

It is not always necessary to displace all the atoms in the primitive cell, since the use of symmetries can reduce the amount of work needed. This is done as follows. We displace one atom in the primitive cell, let's call it 'one', and we calculate the forces induced by the displacement on all the other atoms of the supercell. Then we pick up one other atom of the primitive cell, atom 'two'. If there is a symmetry operation S (not necessarily a point group symmetry operation) such that, when S is applyed to the crystal atom two is sent into atom one and the whole crystal is invariant under such transformation, then it is not necessary to displace atom two, and the part of the force constant matrix associated with its displacement can be calculated using

$${\bf\Phi}_{i s, 0 2} = {\bf B}(S) {\bf\Phi}_{\lambda_{i s}(S), 0 1} {\bf B}(S^{-1}),$$ (8)

where ${\bf B}(S)$ is the $3\times 3$ matrix representing the point group part of S in Cartesian coordinates, and $\lambda_{i s}(S)$ indicates the atom of the crystal where the atom in ${\bf R}_i^0 + {\bf\tau}_s$ is brought because of the action of the symmetry operation S. If there is no symmetry operation connecting atom two to atom one then atom two is displaced and all the induced force field is calculated. The procedure is repeated for all the atoms of the primitive cell.

In principle each atom has to be displaced along the three Cartesian directions. It is sometimes convenient to displace the atoms along some special directions so as to maximize the number of symmetry operations still present in the 'excited' supercell, in this way the calculations of the forces are less expensive. This can always be done, as long as one displaces the atoms along three linearly independent directions. The forces induced by the displacements along the three Cartesian directions is easily reconstructed by the linear combination

$${\bf F_{i s, 0 t \alpha}} = \sum_l {\bf A}_{l \alpha} {\bf \tilde{F}}_{i s,0 t k }$$ (9)

where ${\bf\tilde{F}}_{i s, 0 t k}$ is the force induced on the atom in ${\bf R}_i^0 + {\bf\tau}_s$ due to a displacement of the atom in ${\bf\tau_t}$ along the direction ${\bf u_k}$, and ${\bf A} = {\bf (\frac{u_1} {\mid u_1 \mid},\frac{u_2}{\mid u_2 \mid},\frac{u_3}{\mid u_3 \mid})}^{-1}$is the inverse of the $3\times 3$ matrix whose columns are the normalized displacements in Cartesian coordinates.

Using symmetries it is possible to reduce the number of displacements even further: if applying a point group symmetry operation U to the displacement vector ${\bf u}_1$ one obtains a vector ${\bf u}_2$ which is linearly independent from ${\bf u}_1$, then the force field that would be induced by the displacement ${\bf u}_2$ can be calculated by

$${\bf F}_{i s, 0 t 2 } = {\bf B}(U) {\bf F}_{\lambda_{i s}(U^{-1}), 0 t 1 }.$$ (10)

If a linearly independent direction cannot be found one has to displace the atom along a chosen independent direction and perform an other calculation. This is done untill a set of three independent directions is found.

The force constant matrix is invariant under the point group symmetry operations of the crystal. This is not automatically garanteed by the procedure just described, because in general the crystal is not harmonic, and therefore eqn.  is only an approximation. So, the force constant matrix must be symmetrized with respect to the point group operations of the crystal:

$${\bf\Phi}_{i s, 0 t} = \frac{1}{N_G} \sum_U {\bf B}(U) {\bf\Phi}_{\lambda_{i s}(U), 0 t} {\bf B}(U^{-1}).$$ (11)

The symmetrization of the force constant matrix removes all even-order anharmonicities. The harmonic approximation becomes better and better as the displacement are made smaller and smaller. However, if the displacements are small, also the force induced are small, but there is a limit in the accuracy achievable in the calculations, so one cannot make too small displacements. Usually a fraction of a $\%$ of the nearest-neighbour distance is a good compromise.

As an example of the procedure just described let's consider the h.c.p. crystal. There are two atoms in the primitive cell, so in principles we would need six independent calculations. We will see that the number of calculations needed is equal to two. In first place one can easily recognize that only one atom needs to be displaced: if we traslate the crystal from one atom to the other and we perform a spatial inversion the crystal remains unchanged. Secondly, by applying a clockwise rotation of 120 degrees, for example, to a displacement in the x direction, one obtains an independent displacement. So only one additional displacement along the z direction is needed.