Phonon Calculations for Bilayer 2D systems

Question

I am facing problems in calculating the phonon band structure for bilayer 2D systems. While the band structure mostly shows positive frequencies for the 1 layer system, going to 2 layers makes some frequencies imaginary. Moreover, the choice of van der Waals interaction also comes into account for the bilayer system and for some van der Waals functionals, phonon calculations are not properly implemented.

Can someone share thoughts on calculating the phonon band structure for bilayer or perhaps multilayered 2D systems and which approach should be adopted to tackle this problem?

Answer 1

The presence of imaginary frequencies in a phonon dispersion can have two sources:

1. Physical origin. If the imaginary frequency appears at a $\mathbf{q}$-point in the Brillouin zone that is included in the $\mathbf{q}$-point grid that you explicitly calculate (e.g. one of the points in the $N\times N\times N$ $\mathbf{q}$-point grid if you use a supercell of size $N\times N\times N$), then this means that the imaginary frequency is physical. This is telling you that your structure is not at a local minimum of the potential energy landscape, but instead at a saddle point, and there is a lower energy structure that you should be using instead. So what do you have to do to find this lower energy structure? You have to distort your original structure along the wave vector associated with the imaginary mode to find the real minimum. A previous discussion went into some detail about how to do this.
2. Underconvergence. If the imaginary frequency appears at a $\mathbf{q}$-point in the Brillouin zone that is not included in the $\mathbf{q}$-point grid that you explicitly calculate, then this may be because of underconvergence. When you plot a phonon dispersion you perform a Fourier interpolation over the grid that you explicitly calculate to obtain phonon frequencies at other points in the Brillouin zone that you don't explicitly calculate. This interpolation may introduce spurious imaginary modes if your explicitly calculated $\mathbf{q}$-point grid is not large enough. In this case, the solution is to perform more calculations at larger $\mathbf{q}$-point grids (larger supercells) until either the imaginary frequencies go away, or you can confirm that they are actually physical by explicitly including them in your direct calculation.

So why may you get imaginary frequencies in the bilayer when you didn't get them in the monolayer? For case 1 above, this may be caused by the way in which you constructed the bilayer. If you simply made a copy of the monolayer and placed it next to the original monolayer, then it may be that the interaction between the two layers means that your individual layers want to distort. In this case, the imaginary frequency is telling you which way to distort to get to the energy minimum. For case 2 above, it could be that a converged $\mathbf{q}$-point grid for the monolayer is not converged for the bilayer, which may very well happen in the out-of-plane direction.

Your other question was about using phonons with van der Waals interactions. But this question can be made more general: with what level of theory can phonons be calculated? The answer depends on the method you use to calculate phonons. For finite difference methods: if you can calculate forces, then you can also calculate phonons. As forces are typically available, then you can calculate phonons with most methods using finite differences. If you instead use density functional perturbation theory, here you are limited to the actual implementation that a given code has. For phonons, this is typically restricted to semilocal functionals.