# How to deduce phase transitions from a phonon calculation?

Matter Modeling Asked by koroma on December 9, 2021

I came across the concept of using phonons to establish a material’s dynamical stability, based on whether or not imaginary frequencies are present in its phonon band structure.

What I am struggling with is how to determine the phase transitions based on a phonon band structure with such imaginary frequencies? From what I have read, the recipe is to follow the imaginary phonons. What does this mean? How does one follow such imaginary modes?

An example illustrating this using perovskite transitions from the high symmetry $$Pmoverline{3}m$$ space group to lower symmetries will be appreciated.

The answer of @ProfM is already very complete, but I wanted to tackle your question from a more practical point of view.

The presence of imaginary frequencies indicate that there are atomic positions which are more energetically favorable at the ground state. So, the concept of "following" a mode means condensing it onto the reference structure, until you find the equilibrium positions.

To give an example, we can start with $$ce{BaTiO_3}$$ with $$Pmoverline{3}m$$ symmetry. By looking at the phonon frequencies you will notice the instability (imaginary frequency) at $$Gamma$$, which for this case in particular corresponds to the so-called ferroelectric mode.

Once you know the phonon of interest, you can read its eigendisplacements ($$U_{FE}$$), condense them onto the reference structure ($$S_{ref}$$) using different amplitudes $$alpha$$, $$S_{alpha} = S_{ref} + alpha U_{FE}$$ and then calculate the energy for each of the resulting structures ($$S_{alpha}$$). You will notice that, as you increase the value $$|alpha|$$, the total energy of the system decreases until it starts increasing again, creating a potential well. The structure with the minimum energy will be the one at equilibrium (according to the selected mode).

Coming back to the $$ce{BaTiO_3}$$ example, the new structure $$S_{alpha, E_{min}}$$ should show a $$P4mm$$ symmetry. However, you would need to relax the phase, as strain also needs to be included.

Note: $$alpha$$ can be any real value but it is case-specific (depending on the anharmonicities of the system and on the definition of the eigendisplacements), so would have to try with different options to determine which one is better for your case.

Answered by Alam on December 9, 2021

Background theory. In the harmonic approximation, the potential energy surface (PES) is expanded about an equilibrium point to second order, to obtain the Hamiltonian:

$$hat{H}=sum_{p,alpha}-frac{1}{2m_{alpha}}nabla_{palpha}^2+frac{1}{2}sum_{p,alpha,i}sum_{p^{prime},alpha^{prime},i^{prime}}D_{ialpha;i^{prime}alpha^{prime}}(mathbf{R}_p,mathbf{R}_{p^{prime}})u_{pialpha}u_{p^{prime}i^{prime}alpha^{prime}}.$$

The basic quantity one builds when calculating phonons is the matrix of force constants:

$$D_{ialpha;i^{prime}alpha^{prime}}(mathbf{R}_p,mathbf{R}_{p^{prime}})=frac{partial^2E}{partial u_{pialpha}u_{p^{prime}i^{prime}alpha^{prime}}},$$

which is the expansion coefficient of the second order term in the potential energy surface $$E$$, with $$i$$ labelling the Cartesian direction, $$alpha$$ the atom in the basis, $$mathbf{R}_p$$ is the position of the cell $$p$$ in the crystal, and $$u_{pialpha}$$ is the amplitude of displacement of the corresponding atom. Using periodicity of the crystal, we can define the dynamical matrix at each $$mathbf{q}$$-point of the Brillouin zone as:

$$D_{ialpha;i^{prime}alpha^{prime}}(mathbf{q})=frac{1}{N_psqrt{m_{alpha}m_{alpha^{prime}}}}sum_{mathbf{R}_p,mathbf{R}_{p^{prime}}}D_{ialpha;i^{prime}alpha^{prime}}(mathbf{R}_p,mathbf{R}_{p^{prime}})e^{imathbf{q}cdot(mathbf{R}_p-mathbf{R}_{p^{prime}})},$$

where $$N_p$$ is the number of cells in the supercell over which periodic boundary conditions are applied, and $$m_{alpha}$$ is the mass of atom $$alpha$$. Diagonalizing the dynamical matrix gives eigenvalues $$omega^2_{mathbf{q}nu}$$ and eigenvectors $$v_{mathbf{q}nu;ialpha}$$. From these, it is possible to define a set of normal coordinates:

$$u_{mathbf{q}nu}=frac{1}{sqrt{N_p}}sum_{mathbf{R}_p,i,alpha}sqrt{m_{alpha}}u_{pialpha}e^{-imathbf{q}cdot{mathbf{R}_p}}v_{-mathbf{q}nu;ialpha},$$

in terms of which the Hamiltonian becomes a sum over uncoupled simple harmonic oscillators:

$$hat{H}=sum_{mathbf{q},nu}-frac{1}{2}frac{partial^2}{partial u_{mathbf{q}nu}^2}+frac{1}{2}omega^2_{mathbf{q}nu}u_{mathbf{q}nu}^2.$$

The bosonic quasiparticles labelled by quantum numbers $$(mathbf{q},nu)$$ are called phonons, and have energy $$omega_{mathbf{q}nu}$$ and momentum $$mathbf{q}$$.

Dynamically stable structure. A dynamically stable structure is one whose equilibrium position is at a local minimum of the PES. As such, the eigenvalues $$omega^2_{mathbf{q}nu}$$ of the dynamical matrix (Hessian) are all positive numbers, and as a consequence the phonon frequencies $$omega_{mathbf{q}nu}$$ are all real.

Dynamically unstable structure. A dynamically unstable structure is one whose equilibrium position is at a saddle point of the PES. As such, some of the eigenvalues of the dynamical matrix are negative, and the corresponding phonon frequencies imaginary.

Physical interpretation. Phonons measure the curvature of the PES around the equilibrium position of the material. As we have seen, an imaginary frequency corresponds to a negative curvature, so it corresponds to a direction in the PES along which the energy decreases. This means that there is a lower energy configuration of the material, and we say that the structure is then dynamically unstable.

"Follow the imaginary modes". How can we find such lower-energy structure? The eigenvectors of the dynamical matrix associated with the imaginary phonons tell us the direction along which the energy decreases, so we can "follow those modes" to find the lower energy structure. This can be done by simply constructing a sequence of structures on which you displace the atoms by an amplitude $$u_{mathbf{q}nu}$$ (see equation above) of the imaginary phonon $$(mathbf{q},nu)$$, and calculating the total energy of each of the resulting structures. For a saddle point, the resulting curve will be something like a double well, and the minima of the double well correspond to your new lower-energy structure.

Finite temperature. The discussion up to this point concerns the potential energy surface, so temperature is neglected. If you are interested in a calculation at finite temperature, then you need the free energy surface. This is much harder to calculate and you need anharmonic terms in your Hamiltonian to describe it properly.

Perovskites. Perovskites typically have a cubic structure at high temperature, and then upon lowering the temperature undergo a number of phase transitions to lower-symmetry structures (tetragonal, orthorhombic, and so on). Imagine a perovskite that has only two phases, tetragonal at low temperature and cubic at high temperature (generalizing to more phases is trivial). Then if you calculate the phonons in the cubic structure (saddle point) you will find imaginary modes, and following them will take you to the tetragonal structure (minimum). If you calculate the phonons in the tetragonal structure, they will all have real frequencies. So why is the cubic phase stable at high temperatures? This is because, although the cubic phase corresponds to a saddle point of the potential energy surface, above some critical temperature it corresponds to a minimum of the free energy surface. As such, above that critical temperature the cubic phase becomes dynamically stable. As I mentioned above, to investigate this phase transition (e.g. to calculate the critical temperature), you need to include anharmonic terms (phonon-phonon interactions), which is much harder computationally.

Answered by ProfM on December 9, 2021

## Related Questions

### I get this error in Gaussian09 even though I have enough memory. “Error termination in NtrErr: NtrErr called from FIOCnC.” Any suggestions?

1  Asked on January 5, 2022 by quantumx

### Why are single excitations ignored in the MP2 component of double hybrid functional calculations?

1  Asked on January 5, 2022

### Regarding oscillatory strength theoretical units to experimental ones

1  Asked on January 3, 2022

### How well can we model chemical synthesis?

3  Asked on January 1, 2022 by anyon

### Is it required to sample the full reaction coordinate in umbrella sampling?

1  Asked on January 1, 2022 by bnd

### How to find good universities for matter modeling?

1  Asked on December 25, 2021

### What does the Neighbor command do in LAMMPS?

1  Asked on December 22, 2021

### For software that does not support FCIDUMP format, what format is used and how can we get the software to interact with FCIDUMP integrals?

2  Asked on December 22, 2021

### How do the various programs read or write integrals in FCIDUMP format?

1  Asked on December 20, 2021

### Is there an anisotropy factor (g factor) for TDDFT absorption and circular dichroism calculations?

1  Asked on December 20, 2021 by c-alexander

### Why am I losing bonds when I use PACKMOL on two proteins?

1  Asked on December 20, 2021

### Can we “invert” Density Functional Theory through sufficiently accurate experiment?

1  Asked on December 18, 2021 by kfc

### How do you calculate the “true” chemical potential in classical density functional theory?

1  Asked on December 16, 2021

### Converging the potential of a dilute binary alloy system using the Coherent Potential Approximation. What troubleshooting options could I take?

1  Asked on December 16, 2021 by luphys

### The Hund’s J – Why can this be quantified?

1  Asked on December 16, 2021

### How accurate are the most accurate calculations?

1  Asked on December 13, 2021 by roman-korol

### Is there any software that can generate the electronic DOS of GaAs?

1  Asked on December 13, 2021

### Which Linux distribution is best for Matter Modeling?

9  Asked on December 11, 2021

### What is the role of a chemist in computational chemistry?

4  Asked on December 9, 2021 by georgios-karaiskakis

### How to deduce phase transitions from a phonon calculation?

2  Asked on December 9, 2021 by koroma