How big should a supercell be in phonon calculations?

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.]

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.

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.