4.6 Supercell calculations

The supercell method is the ubiquitous approach for the study of solid-state periodic boundary condition systems. In principle, a solid may also be approximated by means of a large cluster. For sufficiently large clusters, the quantum mechanics of the centre-most atoms approximate those in a solid. However, the size of cluster required to approximate a solid is large due to the dominance of surface over ``bulk'' atoms in small and medium sized clusters. This consideration applies to all electronic structure methods, and the best general solution is to use a limited number of repeats of a primitive cell combined with periodic boundary conditions, as illustrated in figure 4.2.

Figure 4.2: The supercell approach to periodic systems. Periodic boundary conditions are applied and all interactions are periodic with the periodicity of the supercell. The ``minimum image'' convention must be applied when particle-particle distances are computed: the minimum distance between particles can involve particle images in adjacent supercells.

In the supercell approach, [73,17] an artificial periodicity is imposed on the simulation cell to better model the continuum properties of the system. Bloch's theorem may then be applied to the wavefunctions.

For perfect crystals, a limited number of repeats of the primitive cell are used. Periodic boundary conditions impose an artificial periodicity on interactions in the system. The supercell must be large enough that these spurious interactions are small. These finite size effects are examined in detail in chapter 6.

4.6.1 Long-range interactions

Within the supercell approximation, long-ranged interactions require special treatment. Any potential extending significantly beyond a single supercell, such as the coulomb potential, must correctly computed to account for the periodic boundary conditions. Neglect of the contributions to the potential from particles in neighbouring cells can lead to significant errors in energies, and the limiting behaviour of certain properties can be altered. Accurate summation of long-ranged interactions can be computationally intensive. Although approximate schemes have been developed for the rapid evaluation of long-ranged interactions, we are interested in a high final accuracy. Consequently, only accurate schemes for the evaluation of the coulomb interaction have not been used.

Ewald's method [51,74,75,17] is a near exact method for the evaluation of long-ranged interactions within periodic boundary conditions. The interaction terms forming the potential are rearranged into two series that are readily evaluated.

To illustrate the method, we consider a cubic supercell containing $N$ particles of charge $q_{i}$. The interaction energy can be written as

$$E=\frac{1}{2}\sum_{{\bf n}} ^{\prime}\sum_{i=1}^{N}\sum_{j=1}^{N} \frac{q_{i}q_{j}}{\vert{\bf r}_{ij}+{\bf L}\vert} \;\;\;,$$ (4.40)

where ${\bf L}$ are all supercell lattice vectors, and the prime indicates that for ${\bf L}=0$ the self-interaction terms $i=j$ are excluded. This sum is divergent unless the system is overall charge neutral ($\sum q_{i}=0$) [76] and care is required in its evaluation.

In the Ewald method, the point charge distribution is modified to give two convergent series, the sum of which gives the correct value for the unmodified distribution. The neutral Ewald ``charge density'' is given by a sum over the point charges plus a screening background.

$$\rho_{\rm Ewald}({\bf r})=\sum_{{\bf L}}\delta({\bf r}-{\bf r}_i-{\bf L}) -\rho_{\rm background} \;\;\; ,$$ (4.41)

where the ${\bf r}_i$ are the positions the point charges and ${{\bf L}}$ is the complete set of supercell lattice vectors.

The Ewald potential is computed as a sum of two components. An array of Gaussians are used in both: in the first, the Gaussians cancel the background charge density, and in the second, the Gaussians are subtracted from the delta function charge distribution (see figure 4.3). The width of the Gaussian charges, $\mu$, controls the rate of convergence of the two series and is usually optimised for each supercell.

Figure 4.3: The Ewald method for periodic potentials, shown schematically in 1 dimension. A screening charge distribution is added and subtracted to the periodic point charge distribution, giving two rapidly convergent series. Top: The Ewald charge density is a sum of delta functions (red) minus a background charge density (blue). Centre: $\rho_1$ (equation 4.42) is given by an array of Gaussians (magenta) and a background charge density. Bottom: $\rho_2$ (equation 4.43) is given by an array of delta functions minus an array of Gaussians.

$$\rho_1({\bf r}) = \frac{1}{\mu\sqrt{\pi}} \sum_{{\bf L}}\exp({-({\bf r}-{\bf r}_{i}-{\bf L})^2/\mu^2})-\rho_{\rm background}$$ (4.42)

and

$$\rho_2({\bf r}) = \sum_{{\bf L}} \left[ \delta({\bf r}-{\bf ... ...exp({-({\bf r}-{\bf r}_{i}-{\bf L})^2/\mu^2}) \right] \;\;\; ,$$ (4.43)

where $\mu$ is the width of the Gaussian charges. The potential due to $\rho_{1}$ is most easily computed in reciprocal space: a Gaussian in real space becomes a Gaussian in reciprocal space. The potential, $\phi_2$, due to $\rho_{2}$ is given in real space:

$$\phi_2({\bf r}) = \sum_{{\bf L}}\frac{{\rm erfc}(\vert{\bf r... ...f r}-{\bf r}_{i}-{\bf L}\vert} - \frac{\pi\mu^2}{V^2} \;\;\; ,$$ (4.44)

where

$$\mathrm{erfc}(x)=1-\mathrm{erf}(x)= 1-\frac{2}{\sqrt{\pi}}\int_{0}^{x}\exp(-x^{\prime 2})dx^{\prime} \;\;\;$$ (4.45)

and the supercell volume is $V$. The sum runs over supercell lattice vectors, ${\bf L}$, and converges rapidly so that only a few vectors are needed in the sum. The overall Ewald potential is:

$$\phi_{\rm Ewald}({\bf r}) = \sum_{{\bf G}\neq 0} \frac{4\pi}{V{\bf G}^{2}} \exp({-{\mu^2{\bf G}^{2}/4}}+i({\bf r}-{\bf r}_{i}).{\bf G})$$

$$+ \sum_{{\bf L}}\frac{{\rm erfc}(\vert{\bf r}-{\bf r}_{i}-{\bf L}\vert/\mu)} {\vert{\bf r}-{\bf r}_{i}-{\bf L}\vert} - \frac{\pi\mu^{2}}{V} \;\;\;,$$ (4.46)

where ${\bf G}$ are supercell reciprocal lattice vectors.

This potential could be used to compute the electron-ion potential. The Ewald sum for the electron-electron potential is similar, but the ``self-image'' of an electron is removed, due to the periodic array of images and Gaussian charges. Electrostatic energies are computed by multiplying the Ewald potential by the appropriate product of charges $q_{i}q_{j}$. In terms of speed and accuracy, for small numbers of particles (1000 or less), a well implemented version of Ewald's method has yet to be improved on. The computational cost Ewald's method, for a specified accuracy, scales as $N^{3/2}$. [77] Large numbers of particles are more efficiently dealt with using multipole expansions, based on hierarchical spatial subdivision of the charge density. [78,79] The largest supercell calculation in this thesis involves 1000 electrons (see chapter 6), and Ewald's method was applied in this case.