About k-point sampling

Background theory

The fact that the periodicity is an inherent feature of a crystal structure, has very important consequences for the numerical treatment of such structures; these must be considered when we calculate periodic systems using any DFT code.

The fact that a crystal is invariant in certain directions with respect to translations by a unit cell length implies that the wave vector is a conserved quantity (note that, you can have 1D and 2D periodic systems also). As a result, the eigenfunctions obey Bloch theorem and are therefore labeled by the wave vector $\mathbf {k}$ (in 1D, k is just a scalar, the wave number k). This also means that the Schrödinger equation can be solved independently for each possible value of $\mathbf {k}$.

When constructing the total electron density $n(\mathbf {r} )$, it is, however, necessary to know the solution in all points of the first Brillouin zone, since

$n(\mathbf {r} )=\sum _{\alpha }\int d{\mathbf {k} }\,|\Psi _{\alpha \mathbf {k} }(\mathbf {r} )|^{2}f_{n\mathbf {k} },$

where $f_{\alpha \mathbf {k} }$ is the occupation of the state $\mathbf {k}$ in the band $\alpha$ and $\Psi _{n\mathbf {k} }(\mathbf {r} )$ is the corresponding Bloch wavefunction. In practice, the integral is approximated by a sum, which is then evaluated over a finite number of k-point samples in the Brillouin zone. Clearly, as the number of k-points increases, the sum approximates the "true" value of the above integral to a greater accuracy.

Exactly how many k-points to choose depends to a large extent on the geometry and dimensionality of the system. For example, the bigger the unit cell is, the fewer points are generally needed. Please note, that it is only necessary to sample the directions in which the crystal is periodic. A carbon nanotube, for example, only requires a large k-point sampling along the tube axis (and a single point, of course, in the two transverse directions).

Functions for k-point sampling in GPAW

In GPAW, the default sampling of the Brillouin zone is with only the $\Gamma$-point. This allows us to choose the wave functions to be real. Monkhorst-Pack sampling can be used if required:

kpts=(N1, N2, N3)

where N1, N2 and N3 are positive integers. This will sample the Brillouin zone with a regular grid of $N1\times N2\times N3$ k-points. See the ase.dft.kpoints.monkhorst_pack() function for more details.

For more flexibility, you can use this syntax:

kpts={'size': (4, 4, 4)}  # 4x4x4 Monkhorst-pack
kpts={'size': (4, 4, 4), 'gamma': True}  # shifted 4x4x4 Monkhorst-pack

You can also specify the k-point density in units of points per Ang^-1:

kpts={'density': 2.5}  # Monkhorst-Pack with a density of 2.5 points/Ang^-1
kpts={'density': 2.5, 'even': True}  # round off to neares even number
kpts={'density': 2.5, 'gamma': True}  # include gamma-point

The k-point density is calculated as $Na/2\pi$, where $N$ is then number of k-points and $a$ is the length of the unit-cell along the direction of the corresponding reciprocal lattice vector.

An arbitrary set of k-points can be specified, by giving a sequence of k-point coordinates like this:

kpts=[(0, 0, -0.25), (0, 0, 0), (0, 0, 0.25), (0, 0, 0.5)]

The k-point coordinates are given in scaled coordinates, relative to the basis vectors of the reciprocal unit cell. The above four k-points are equivalent to:

kpts={'size': (1, 1, 4), 'gamma': True} 

and to this:

from ase.dft.kpoints import monkhorst_pack
kpts = monkhorst_pack((1, 1, 4))
kpts
array([[ 0.   ,  0.   , -0.375],
      [ 0.   ,  0.   , -0.125],
      [ 0.   ,  0.   ,  0.125],
      [ 0.   ,  0.   ,  0.375]])
kpts+=(0,0,0.125)
kpts
array([[ 0.  ,  0.  , -0.25],
      [ 0.  ,  0.  ,  0.  ],
      [ 0.  ,  0.  ,  0.25],
      [ 0.  ,  0.  ,  0.5 ]])