Finite displacement methods

The capability to obtain forces on ions allows simple calculation of phonons at specific parts of the Brillouin Zone, according to the number of primitive unit cells included in the supercell. Phonons are calculated from the dynamical matrix, which contains spring constants ${\bf\phi}$ representing the interactions between pairs of atoms. An atom $\kappa$ in primitive cell l, having mass $m_\kappa$, is displaced from equilibrium by ${\bf u}\left(^{l}_{\kappa}\right)$ as the lattice vibrates. The force on each atom is then given by

$$F_{\alpha} \left(^{l}_{\kappa}\right) = \sum_{l'} \sum_{\kap... ...{l'}_{\kappa'} \right) u_{\beta} \left(^{l'}_{\kappa'}\right)$$ (4.1)

or, alternatively, as the second derivative of the total energy E:

$$\Phi_{\alpha \beta}\left( ^{l}_{\kappa} \ ^{l'}_{\kappa'} \ri... ...\right) \, \partial u_{\beta}\!\left(^{l'}_{\kappa'}\right)}.$$ (4.2)

A normal mode with angular frequency $\omega$ is defined as a motion such that ${\bf F}(l,\kappa)= - m_{\kappa} \omega^2 {\bf u}(l,\kappa)$. We naturally expect phonon solutions, in which two equivalent atoms, having the same $\kappa$ but in different primitive cells, will differ in the phase of their motion by an amount ${\bf q}\!\cdot\!\Delta{\bf x}$, where $\Delta{\bf x}$ is the difference in the positions of their primitive cells. We define mass reduced coordinates $\mbox{\boldmath$\epsilon$ } \left(^{l}_{\kappa}\right) = \sqrt{m_{\kappa}} {\bf u}\!\left(^{l}_{\kappa}\right)$ and, using (4.1), the equations of motion for phonon solutions define an eigenproblem at each wavevector q:

$$- \omega^2 \mbox{\boldmath$\xi$ } = {\bf D}(\bf q) \mbox{\boldmath$\xi$ }$$ (4.3)

where the Fourier-transformed dynamical matrix is

$$D_{\alpha \beta}(\kappa \kappa'\vert{\bf q}) = \sum_{l} \fra... ...^{l}_{\kappa'} \right) \exp(i{\bf q}.[{\bf x}_l-{\bf x}_0]) .$$ (4.4)

For each wavevector q, evaluation and diagonalisation of ${\bf D(q)}$will give the frequencies and eigenvectors $\xi$ of normal modes at q, which may be assigned to the corresponding point of the Brillouin Zone.

In principle, interactions between pairs of atoms at all separations up to the infinite are required, but in practice the non-Coulombic contribution decays rapidly with distance. Furthermore, it can be shown that in a supercell with periodic boundary conditions, calculations of a dynamical matrix for all atom pairs within one supercell is sufficient to give phonons at all wavevectors commensurate with the supercell, as forces are calculated to include all long-range contributions over other supercells. For example, calculation of phonons at $\Gamma$ may be done with only a single primitive cell.

The dynamical matrix for the supercell may be found by making small displacements of one atom at a time, and calculating the forces exerted on all atoms. From (4.1) the force constants may thus be calculated. In a system with symmetry, the symmetry elements may used to deduce related force constants, thus minimising the number of ab initio simulations that need to be performed (14).

An example is shown in Table 4.1, for calculations on two phases of SiO₂, related by a phase transition at 47 GPa (from (15)). For many modes the agreement between theory and experiment is better than a few cm $^{\scriptscriptstyle -1}$. Four different displacements were needed for the rutile phase, but six for the lower-symmetry CaCl₂ phase; for each an average force constant over positive and negative displacements was taken.

 

Table: LDA frequencies (cm $^{\scriptscriptstyle -1}$) of Raman TO modes at $\Gamma$ in SiO₂ in the stishovite (rutile structure) phase at 0 GPa, and CaCl₂ phase at 50 GPa, with finite-displacement (DFT) [15] and linear-response (LR) [23] methods, compared to available experimental data (references in [15]). Similar agreement is found for infrared modes.

Rutile				CaCl2
	DFT	LR	Expt.		DFT	Expt.
Raman
B1g	222	214	232	Ag	174	163a
Eg	582	585	589	B2g	669	684
				B3g	675	706
A2g	598	599	Silent	B1g	673
A1g	750	755	754	Ag	889	906
B2g	947	954	966	B1g	1111
a Broadband frequency is 177 cm $^{\scriptscriptstyle -1}$
b Transmission and reflection spectra differ

 

Diamond-structured silicon has extremely high symmetry, with a primitive unit cell of only two atoms, so it is straightforward to construct a supercell containing many primitive unit cells and thus obtain phonons at several points in the Brillouin Zone. Force constants were obtained for a cubic cell containing 64 atoms, from a single calculation in which one atom was slightly displaced. In this and other systems it has been found that enforcing the reduced symmetry of the distorted structure significantly reduces the noise in the calculations (16), although anharmonicity may still be present.

Since the interactions in Si decay rapidly with distance, to a good approximation they may be neglected for a distance of more than half the supercell, and phonons at any wavevector may then be calculated. The dispersion curve in Figure 4.1 was obtained from these calculations, which agrees well with experiment. The small amount of degeneracy breaking at W is due to symmetry breaking in the original calculated force constants, but is hard to enforce without compromising results elsewhere in the Brillouin Zone (16).

  

Figure 4.1: Dispersion curve of Si calculated from displacement of one atom in a 64-atom supercell. Lines are LDA calculations [16]; points are from neutron scattering experiments (references in [17]).

Apart from the restriction on wavevectors, the other major limitation of this method is that periodic boundary conditions do not permit a dipole in the unit cell, as this would necessitate a non-periodic potential. At the zone centre, longitudinal optic (LO) phonons in polar materials produce such a dipole, and thus are not correctly calculated; they are instead found to be degenerate with transverse optic (TO) modes.

LO and TO modes at $\Gamma$ are related by the Lyddane-Sachs-Teller relation, which for a structure with only two atoms in the primitive cell takes the form

$$\frac{\omega_{LO}^2}{\omega_{TO}^2} = \frac{\epsilon_0}{\epsilon_\infty} .$$ (4.5)

For such a structure, the values of the dielectric constants are sufficient to obtain the corrected LO frequency. In more complex structures they are included in the Fourier-transformed dynamical matrix with additional terms taking the form (17)

$$D'_{\alpha \beta}(\kappa \kappa'\vert{\bf q}) = \frac{4 \pi e... ...bf q} \cdot \mbox{\boldmath$\epsilon$ }^\infty \cdot {\bf q} }$$ (4.6)

so the high-frequency static dielectric tensor $\epsilon$$^\infty$ and the Born effective charge tensors ${\bf Z}^\ast_\kappa$ are needed. These are not straightforward to calculate, but are more naturally found using the formalism outlined in the next section.