Electro-elastic model for phonon calculations

The vibrations of atoms in a solid material around their equilibrium positions give rise to the quantum-mechanical phonons of the system. Unfortunately, phonon calculations in the DFT framework are very costly, in particular when phonons must be calculated throughout the Brillouin zone of a crystal. It is therefore common to set up an empirical model for the interatomic interactions and calculate the phonons from the model, and use DFT calculations to determine the parameters in the model. In this project, we introduced a mixed Coulomb-Hook model with anisotropic Born effective charges in order to calculate phonons in ionic systems with defects.

Motivation

MgO phonon bandstructure

The vibrations of atoms in a solid material around their equilibrium positions give rise to the quantum-mechanical phonons of the system. Fortunately, the vibrational frequency in the quantum-mechanical and classical theory are the same in the harmonic approximation, i.e., when only the curvature of the energy surface at the equilibrium positions are considered. This gives a direct link between the phonon modes in the system and the total energy surface as calculated e.g. from density-functional theory, our work horse in ab initio electronic structure theory. Yet, phonon calculations in the DFT framework are very costly, in particular when phonons must be calculated throughout the Brillouin zone of a crystal.

Coulomb-Hook model

Electro-elastic interaction model. Each ion bears a charge and interacts with all other ions via the Coulomb interaction. In addition, short-range chemical interactions (depicted by springs) occur between neighbouring atoms.

However, a full DFT treatment of the phonons is rarely needed. Since the phonons only depend on the interatomic interactions, and these rapidly decay with interatomic distance, it is sufficient to set up an empirical model for the interatomic interactions and calculate the phonons from the model, and use DFT calculations to determine the parameters in the model. The most simple approach is the truncated Hook model, where we introduce a spring constant between nearby atoms. The number of nearby atoms to be included depends - of course - on the decay of the interactions with distance. The model therefore works nicely for metals, but fails in practice for ionic materials where Coulomb interactions between the charged ions mediate long-range interactions. The basic idea of the Coulomb-Hook model, therefore, is to use charge-charge Coulomb interactions in addition to the truncated Hook model. The idea is not new, and has been used in the past to calculate phonons in bulk materials. It is also formally well founded, since the exact DFT phonons can be shown to assume the analytic form of the Coulomb-Hook model. The charges of the model are known as Born effective charges, and can be defined as the change in the system dipole when an atom is displaced from its equilibrium position.

Anisotropic Born charges

What is new in our approach is that we consider anisotropic Born effective charges. They arise when the apparent charge depends on the direction in which the atom is moved. This regularly happens near point defects, where movement towards the defect center can be very different from other directions. The generalization of an isotropic charge model to anisotropic charges is far from trivial when we require that important fundamental properties should be maintained: notably the existence of a well defined energy surface (which implies the symmetry of the Hesse matrix) and translational symmetry (when all atoms move together, the energy must not change).

To solve this issues, we derive our Coulomb interactions from charges and additional dipoles centered on the atoms. The dipoles, in turn, depend on the interatomic distances. We have shown that such a model fulfills the required symmetry properties and is also general enough to allow for a one-to-one mapping to DFT. A systematic way to achieve the necessary parametrization from a few DFT force calculations has also been devised. Last, but not least, we have gone through the lengthy derivation of the analytical Hesse matrix and implemented an efficient calculation using Ewald summation.

The code has been used to calculate phonons of the oxygen vacancy in MgO.

Go to Editor View