# Are there differences in accuracy and reliability between the frozen phonon method and Density Functional Perturbation Theory?

When proposing or predicting new materials, whether in the bulk or in two dimensions, it is important to discuss the stability of the system. This is usually done through thermodynamics and(or) lattice dynamics. A feasible way to determine the dynamical stability of a material is by modeling its phonon band structure.

I see that two popular methods of modeling the phonon band structure are:

(a) the frozen phonon / finite displacement method (as implemented in phonopy)

(b) Linear Response Theory / Density Functional Perturbation Theory (as implemented in Quantum ESPRESSO (ph.x)

Are the two methods equivalent or is one more reliable/accurate than the other, and if so, why?

Short answer: Modern implementations of these two methods lead to similar accuracies.

Longer answer: The calculation of phonons requires the calculation of the Hessian of the potential energy surface $$V(\mathbf{R})$$, also known as the matrix of force constants:

$$\frac{\partial^2 V(\mathbf{R})}{\partial \mathbf{R}_i\partial\mathbf{R}_j}=-\frac{\partial \mathbf{F}_j}{\partial\mathbf{R}_i},$$

where $$\mathbf{R}$$ is a collective coordinate of all atomic positions, $$i$$ and $$j$$ label atoms $$i$$ and $$j$$ in your system, and $$\mathbf{F}_j=-\partial V(\mathbf{R})/\partial\mathbf{R}_j$$ is the force felt when displacing atom $$j$$.

Finite displacement. This is what you call "frozen phonon", and in this method you calculate the forces in DFT, and then calculate the derivative of the forces by finite difference methods. Therefore, the numerical approximation is that of approximating a derivative with a finite difference formula. In principle, you can make this calculation as accurate as you want by using increasingly accurate approximations to the numerical derivative, but in practice even low-order approximations lead to very accurate answers. The advantages of this method are that it is very simple to implement, and therefore it is in fact available using any underlying electronic structure method that can calculate forces, which includes semilocal DFT, hybrid DFT, or other non-DFT methods, like force fields or dynamical mean-field theory. The disadvantage of this method is that it requires the construction of supercells to capture long-wavelength phonons, which can make the calculations expensive. Most finite displacement codes use "diagonal" supercells, which lead to poor scaling, but the recently introduced "nondiagonal" supercells here [disclaimer: I am a co-author of this work], significantly reduce the computational cost.

DFPT. In density functional perturbation theory, the calculation of the Hessian is specialized to DFT, and the second derivative of the energy is calculated as

$$\frac{\partial^2 E}{\partial\lambda_i\lambda_j}=\int\frac{\partial^2V(\mathbf{r})}{\partial\lambda_i\lambda_j}n(\mathbf{r})d\mathbf{r}+\int\frac{\partial n(\mathbf{r})}{\partial\lambda_i}\frac{V(\mathbf{r})}{\partial\lambda_j}d\mathbf{r}.$$

This expression is general for parameters $$\lambda$$, and in the case of phonons they simply become the atomic coordinate $$\lambda_i=\mathbf{R}_i$$. This expression requires the calculation of the derivative of the density, which in turn requires the derivative of the Kohn-Sham states. These terms can all be calculated within DFPT with the usual numerical approximations of finite basis sets, etc. Again, in practice, modern implementations are relatively easy to converge. The advantage of DFPT is that it does not require the construction of supercells, one can build a finite wave vector response within the primitive cell, so the computational cost is smaller than in the finite displacement method. The disadvantage of this method is that it is restricted to DFT (so no DMFT for example), and furthermore, the algorithmic implementation is not trivial, so it is only widely available for semilocal DFT (so no hybrid DFT either).

In summary, these two methods lead to comparable accuracies. If DFPT is available, then the calculations will be cheaper and DFPT should be the method of choice. However, DFPT is only widely available with semilocal DFT, which means that if you want phonons at the hybrid functional level, or using beyond-DFT methods, then you have to use the finite displacement method.

• Just a small comment: DFPT for hybrid functionals is actually very widely available in molecular codes (some codes support even double hybrid functionals), but in molecular codes it is usually called an analytic Hessian (or analytic frequency) calculation instead of DFPT, and the results are called vibrational frequencies and normal modes, instead of phonons. Therefore, if the desired phonon modes are well described by a cluster model, then one can do hybrid functional-based DFPT calculations with a molecular code. Apr 10, 2022 at 12:42