# How big should a supercell be in phonon calculations?

Matter Modeling Asked on August 19, 2021

Is 2x2x2 enough, or is larger needed? I know that a convergence test must be carried out, but increasing supercell size enormously increases computational time, and I am using the Python program Phonopy.

Considering this in terms of a 2x2x2 supercell matrix is the wrong way to think about this, as choice depends on the cell length and bonding type. Given that rigorous convergence testing is near-impossible (see ProfM's answer), what saves the method is the rapid fall-off with distance of the force constant matrix $$Phi$$. This "nearsigtedness" means that the effect of a displacement of one atom will decrease with distance with a rapid power-law (see Gonze et al). They show the decay of force constants in Quartz to near-zero at a distance of ~$$7.5 overset{circ}{mathrm{A}}$$. Consequently a cutoff radius of $$7.5 overset{circ}{mathrm{A}}$$, or a diameter of $$15 overset{circ}{mathrm{A}}$$ will suffice.

So in that case a supercell dimension of $$15 overset{circ}{mathrm{A}}$$ on a side or greater will be large enough to ensure that there is no aliasing error from overlap of the (supercell)-periodic images of the force constant matrix [*]. The 7A referred to in a previous answer is unlikely to be sufficient. This is strongly material/bonding-dependent. In polar materials and very rigid materials the fall-off is slower than in soft materials. For example in graphite the in-plane force constants decay much more slowly than the interlayer ones.

[*] Gonze et al. also show that a model correction for the Coulomb dipole contribution can be used to reduce the acceptable cutoff, but this is not usually done for finite-displacement calculations as Born charges and dielectric permittivity tensors are required.]

Answered by Keith Refson on August 19, 2021

Quick Summary: There's no way around performing a convergence test. However, it is possible to obtain convergence much faster than the Phonopy approach by using nondiagonal supercells [1].

The basic quantity you build when performing a phonon calculation is the matrix of force constants, given by:

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

where $$E$$ is the potential energy surface in which the nuclei move, $$u_{palpha i}$$ is the displacement of atom $$alpha$$ (of all atoms in the basis), in Cartesian direction $$i$$ ($$x$$, $$y$$, $$z$$), and located in the cell within the supercell at $$mathbf{R}_p$$. This matrix of force constants is, roughly speaking, measuring the following: if I move an atom at $$mathbf{R}_p$$, what force does an atom at $$mathbf{R}_{p^{prime}}$$ feel? If the atoms are sufficiently far apart $$|mathbf{R}_p-mathbf{R}_{p^{prime}}|gg1$$, then the atoms do not feel the force, and $$D_{ialpha,i^{prime}alpha^{prime}}(mathbf{R}_p,mathbf{R}_{p^{prime}})to0$$. So you need a supercell large enough to capture all relevant non-zero entries in the matrix of force constants. An equivalent picture emerges when we consider the relation between a supercell of size $$N_1times N_2times N_3$$, which is equivalent to sampling the Brillouin zone (BZ) of the system with a $$mathbf{q}$$ grid of size $$N_1times N_2times N_3$$ (as the BZ of the supercell is correspondingly smaller compared to the BZ of the primitive cell). In this language, you need a $$mathbf{q}$$-point grid sampling the BZ that is large enough.

So how quickly does $$D_{ialpha,i^{prime}alpha^{prime}}(mathbf{R}_p,mathbf{R}_{p^{prime}})$$ go to zero? There is no general answer to this question, it is system-dependent. Therefore, you must perform a convergence test. One thing to note is that the size of your primitive cell will play a role: if you are looking at diamond, with a very small primitive cell containing only 2 atoms, then a $$2times 2times 2$$ grid will definitely not be large enough. However, if you consider a system with a primitive cell containing many atoms, for example $$ce{In_2O_3}$$ with 40 atoms in the primitive cell, then a $$2times 2times 2$$ grid may be enough. Another thing to consider is the shape of the primitive cell. If your primitive cell is very elongated along one direction, then the distances are larger already along that direction, so you would probably be better off with a non-uniform sampling grid.

Diagonal supercells. So how are calculations performed in practice? When you need to sample a $$mathbf{q}$$-point grid of size $$N_1times N_2times N_3$$, then a code like Phonopy builds a supercells of size $$N_1times N_2times N_3$$. This is accomplished using what I call a diagonal supercell:

$$begin{pmatrix} mathbf{a}_{s_1} \ mathbf{a}_{s_2} \ mathbf{a}_{s_3} end{pmatrix} = begin{pmatrix} N_1 & 0 & 0 \ 0 & N_2 & 0 \ 0 & 0 & N_3 end{pmatrix} begin{pmatrix} mathbf{a}_{p_1} \ mathbf{a}_{p_2} \ mathbf{a}_{p_3} end{pmatrix},$$ where $$(mathbf{a}_{s_1},mathbf{a}_{s_2},mathbf{a}_{s_3})$$ are the supercell lattice parameters, and $$(mathbf{a}_{p_1},mathbf{a}_{p_2},mathbf{a}_{p_3})$$ are the primitive cell lattice parameters. As you correctly say, this can very quickly become extremely expensive computationally. There are many published phonon calculations that are not properly converged because of this computational bottleneck. However, if you want to do a proper job, there is no way around performing a convergence test. However, things can be done better than this.

Nondiagonal supercells. It was recently pointed out that, in order to sample a $$mathbf{q}$$-point grid of size $$N_1times N_2times N_3$$, it is possible to build smaller supercells that are mathematicall exactly equivalent to the diagonal supercells. These are called nondiagonal supercells because they exploit the fact that you can build an equally valid supercell not only by scaling the primitive cell lattice vectors, but also by making linear combinations of them. In this case, you obtain:

$$begin{pmatrix} mathbf{a}_{s_1} \ mathbf{a}_{s_2} \ mathbf{a}_{s_3} end{pmatrix} = begin{pmatrix} S_{11} & S_{12} & S_{13} \ S_{21} & S_{22} & S_{23} \ S_{31} & S_{32} & S_{33} end{pmatrix} begin{pmatrix} mathbf{a}_{p_1} \ mathbf{a}_{p_2} \ mathbf{a}_{p_3} end{pmatrix},$$

where the $$S_{ij}$$ entries are not necessarily zero for $$ineq j$$. Exploiting this additional degree of freedom, then when you want to sample a $$mathbf{q}$$-point grid of size $$N_1times N_2times N_3$$, the largest supercell you need is of size given by the lowest common multiple of $$N_1$$, $$N_2$$, and $$N_3$$.

This leads to a dramatic reduction of computational time: if you are interested in sampling a $$mathbf{q}$$-point grid of size $$Ntimes Ntimes N$$, then with diagonal supercells (e.g. Phonopy), you need a supercell of size $$N^3$$. With nondiagonal supercells, you need a supercell of size $$N$$. In the original paper there is an extreme example for calculating the phonons of diamond using a $$mathbf{q}$$-point grid of size $$48times48times48$$. Using Phonopy this would be completely impossible, as it would require a supercell of size 110,592 (containing 221,184 atoms)! This calculation is in fact possible (and relatively easy), using nondiagonal supercells, which would only require a supercell of size 48 (containing 96 atoms).

Disclaimer: I am an author of the nondiagonal supercell paper.

1. Lloyd-Williams, J., & Monserrat, B. (2015). Lattice dynamics and electron-phonon coupling calculations using nondiagonal supercells, Phys. Rev. B, 92, 184301 DOI: 10.1103/PhysRevB.92.184301.

Answered by ProfM on August 19, 2021

Ideally, a convergence test would be the best way to decide the required size of the supercell, but it can get expensive.

When phonopy (or any similar calculation technique) finds displacements in the cell based on symmetry, the idea is to see how the displacement of certain ions affects the forces on every ion within the cell. We must then take care that the displaced ions don't affect themselves i.e. the effect of their displacement doesn't interact with itself. I think every type of atom would exert a different potential, which might also depend on the # of electrons taken to be valence electrons, so there can't be a definitive answer to this question (except convergence).

Luckily, on the Sourceforge forum for phonopy they suggest using a supercell with at least 7 angstrom length in each direction.

Answered by Hitanshu Sachania on August 19, 2021

## Related Questions

### What are typical RPA capabilities of plane-wave codes?

1  Asked on August 19, 2021 by michael-f-herbst

### Quadrature over three Euler Angles for orientation averaging

1  Asked on August 19, 2021 by emil-zak

### What are some open-source all-electron DFT alternatives to Wien2K?

4  Asked on August 19, 2021

### Real space projection vs reciprocal space projection in DFT calculations

1  Asked on August 19, 2021

### How to calculate the Fock matrix in the molecular orbital basis PySCF?

1  Asked on August 19, 2021 by wychh

### Bond Order: When and how is it used today?

2  Asked on August 19, 2021

### Is it possible to calculate/estimate the value of the J parameter to be used in the Heisenberg/Ising Hamiltonians?

1  Asked on August 19, 2021

### What are some methods for modeling bulk phase infrared spectra?

2  Asked on August 19, 2021

### What is the closest thing we have to “the” universal density functional?

4  Asked on August 19, 2021

### Molden installation on macOS?

1  Asked on August 19, 2021 by dmitry-eremin

### What are the types of bond orders?

4  Asked on August 19, 2021

### What is band inversion and how to recognize it in band structure?

1  Asked on August 19, 2021

### What are the types of charge analysis?

10  Asked on August 19, 2021 by dgkang

### What are the types of DFT?

5  Asked on August 19, 2021

### Supercomputers around the world

14  Asked on August 19, 2021

### What are the types of pseudopotentials?

2  Asked on August 19, 2021

### What are the types of ab initio Molecular Dynamics?

3  Asked on August 19, 2021 by etienne-palos

### Are there any MD packages that do proper free energy sampling in the NPT ensemble?

1  Asked on August 19, 2021 by b-kelly

### VASP vs MedeA-VASP

1  Asked on August 19, 2021

### Derivation of Slater-Koster equations

2  Asked on August 19, 2021

### Ask a Question

Get help from others!