Choice of the k-point mesh

  For a periodic system, integrals in real space over the (infinitely extended) system are replaced by integrals over the (finite) first Brillouin zone in reciprocal space, by virtue of Bloch's theorem. In fhi98md, such integrals are performed by summing the function values of the integrand (for instance: the charge density) at a finite number of points in the Brillouin zone, called the k-point mesh.  Choosing a sufficiently dense mesh of integration points is crucial for the convergence of the results, and is therefore one of the major objectives when performing convergence tests. Here it should be noted that there is no variational principle governing the convergence with respect to the k-point mesh. This means that the total energy does not necessarily show a monotonous behavior when the density of the k-point mesh is increased.

 

• Monkhorst-Pack mesh
  In order to facilitate the choice of k-points, the fhi98md package offers the possibility to choose k-points according to the scheme proposed by Monkhorst and Pack [30]. This essentially means that the sampling k-points are distributed homogeneously in the Brillouin zone, with rows or columns of k-points running parallel to the reciprocal lattice vectors that span the Brillouin zone. This option is enabled by setting t_kpoint_rel to .true., which should be the default for total energy calculations. The density of k-points can be chosen by the folding parameters i_facs(1..3). With these parameters, you specify to cover the entire Brillouin zone by a mesh of $i\_facs(1) \times i\_facs(2) \times i\_facs(3) \times nkpt$ points. The details of this procedure are as follows: In fhi98md, the Brillouin zone is spanned by the reciprocal lattice vectors ${\bf b}_1,\;{\bf b}_2$ and ${\bf b}_3$ attached to the origin of the coordinate system. According to this definition, one corner of the Brillouin zone rests in the origin. The entire Brillouin zone is tiled by small polyhedra of the same shape as the Brillouin zone itself. The parameters specify how many tiles you have along the ${\bf b}_1,\;{\bf b}_2$ and ${\bf b}_3$ direction. In each tile, you specify k-points supplied in form of a list. The coordinates of these k-points are given relative to the spanning vectors of a small polyhedron or 'tile', i.e.
  
  $${\bf k} = xk(1) {\bf b}_1 + xk(2) {\bf b}_2 + xk(3) {\bf b}_3$$
  
  
  The supplied k-point pattern is then spread out over the whole Brillouin zone by translations of the tile. In other words, the k-point pattern of a smaller Brillouin zone (which would correspond to a larger unit cell in real space) is 'unfolded' in the Brillouin zone of your system under study. Normally, the pattern consists only of a single point in the center of the tile, leading to the conventional Monkhorst-Pack k-point sets.
  • k-point set for a bulk calculation
    A k-point set typically used in a bulk calculation could look like

  Parameter	Value
  nkpt	1	number of k-points supplied
  xk(1..3),wkpt	0.5 0.5 0.5 1.0	k-points and weights
  i_facs(1..3)	4 4 4	k-point folding factors
  t_kpoint_rel	.true.	frame of reference for k-points xk(1..3)

  • k-point set for a slab calculation
    For a surface calculation with the z-axis as the surface normal, you want the k-point mesh to lie in the xy-plane. There is no dispersion of the electronic band structure of the slab in z-direction to sample. If there would be, it just means that the repeated slabs are not decoupled as they should be, i.e. the vacuum region was chosen too thin. Therefore the z-coordinate of all k-points should be zero. The input typically looks like

    Parameter	Value
    nkpt	1	number of k-points supplied
    xk(1..3)	0.5 0.5 0.0 1.0	k-points and weights
    i_facs(1..3)	8 8 1	k-point folding factors
    t_kpoint_rel	.true.	frame of reference
    		for k-points xk(1..3)

  
  

   

  Figure: 2D Brillouin zone of a surface with cubic symmetry with a $8\times 8$ Monkhorst-Pack grid. The thin square indicates the conventional first Brillouin zone, the thick square marks the Brillouin zone as realized in the fhi98md code. The location of one special k-point (out of 64) within its tile is marked by the cross.

  
  
  

  Note: We recommend to use even numbers for the folding parameters. As a general rule, one should avoid using high symmetry points in the Brillouin zone as sampling points, because this would result in an inferior sampling quality at comparable numerical effort, compared to a similar number of off-axis k-points. Conventionally (in contrast to our above definition), the Brillouin zone is chosen to have the origin in its center. For odd numbers of the folding parameters and the setting '0.5 0.5 ...', some of the 'unfolded' k-points will fall on the zone boundary of the conventional Brillouin zone, which is often a symmetry plane. Likewise, the k-point set may contain a periodic image of the $\Gamma$-point. This is normally undesirable.

• The concept of equivalent k-points
  Usually one is not interested in the total energies themselves, but in comparing different structures, i.e. accurate energy differences are required. If the two structures have the same unit cell, the comparison should always be done using the same k-point set, so that possible errors from a non-converged k-point sampling tend to cancel out. A similar strategy can also be applied when comparing structures with different unit cells. We allude to this concept here as 'equivalent k-point sampling': The structure with a large unit cell has a smaller Brillouin zone associated with it. The k-points sampling along this smaller Brillouin zone should be chosen as a subset of the k-point mesh in the larger Brillouin zone, such that the position of the k-points in this subset, expressed in Cartesian coordinates in reciprocal space, agree in both calculations (to check whether this is actually the case, inspect the list of k-points appearing in the inp.ini file). This goal can be achieved in a simple way by choosing appropriate i_facs. As an example, imagine you want to compare two slab calculations, one with a $(2 \times 1)$, the other with a $(4 \times 2)$ unit cell. In this case, use

  Parameter	Value
  i_facs(1..3)	4 8 1	k-point folding factors

  in the first case, and

  Parameter	Value
  i_facs(1..3)	2 4 1	k-point folding factors

  in the second case, leaving the other parameters unchanged.

  Note: ${\bf b}_1$ is orthogonal to the real lattice vectors ${\bf a}_2$ and ${\bf a}_3$. If ${\bf a}_1$ is the long edge of your real space unit cell, ${\bf b}_1$ spans the short edge of your Brillouin zone. Therefore, the k-point sampling mesh has fewer points in the ${\bf b}_1$ direction and more points in the ${\bf b}_2$ direction in the above example.

   

• Chadi-Cohen mesh
  Another convention for choosing a k-point mesh has been proposed by Chadi and Cohen [31], and has been applied to slab calculations by Cunningham[32]. In contrast to Monkhorst and Pack, the refinement of the k-point mesh to obtain higher sampling density is based on a recursive scheme. However, for cubic symmetry, the outcome of this algorithm can also be interpreted as a special Monkhorst-Pack grid. To discuss differences between the schemes, we resort to the simple case of a two-dimensional mesh for a slab calculation. An example where Cunningham's scheme leads to results different from Monkhorst-Pack are systems with hexagonal symmetry, e.g. slabs with (111) surface of fcc-metals. Here, Cunningham proposes to use a hexagonal k-point mesh. To realize such meshes in the fhi98md code, one has to provide explicitly a list of k-points forming the desired pattern. Cunningham's 6-point pattern in the full Brillouin zone can be obtained as follows

  Parameter	Value
  nkpt	6	number of k-points
  xk(1..3),wkpt	0.33333 0.00000 0.0 0.16667	k-points
  xk(1..3),wkpt	0.66667 0.00000 0.0 0.16667	and weights
  xk(1..3),wkpt	0.00000 0.33333 0.0 0.16667
  xk(1..3),wkpt	0.00000 0.66667 0.0 0.16667
  xk(1..3),wkpt	0.66667 0.33333 0.0 0.16667
  xk(1..3),wkpt	0.33333 0.66667 0.0 0.16667
  i_facs(1..3)	1 1 1	k-point folding factors
  t_kpoint_rel	.true.	frame of reference for k-points

  
  

   

  Figure 3.2: 2D Brillouin zone of a fcc(111) surface with hexagonal symmetry with set of 6 special k-points following Cunningham. The thin polygon indicates the conventional first Brillouin zone, the thick polygon marks the Brillouin zone as realized in the fhi98md code.

  
  
  

  When a denser mesh in the same cell is desired, Cunningham's 18-point pattern is obtained from the input

  Parameter	Value
  nkpt	2	number of k-points supplied
  xk(1..3),wkpt	0.3333 0.3333 0.0 0.5	k-points and weights
  xk(1..3),wkpt	0.6666 0.6666 0.0 0.5	k-points and weights
  i_fa(1..3)	3 3 1	folding factors
  t_kpoint_rel	.true.	frame of reference
  		for k-points

  
  

   

  Figure 3.3: 2D Brillouin zone of a fcc(111) surface with hexagonal symmetry with set of 18 special k-points following Cunningham. The thin polygon indicates the conventional first Brillouin zone, the thick polygon marks the Brillouin zone as realized in the fhi98md code.

  
  
  

  Here we have made use of the 'tiling' strategy  employed in fhi98md. An even denser k-point set, consisting of 54 points in the full Brillouin zone, may be obtained by using the 6 k-points of the first example, but as a pattern repeated in each of the nine tiles, i.e. by setting the folding parameters in the first example to 3 3 1.

• User-supplied k-point sets
  In some cases (like a band structure calculation), the user might find it more convenient to specify the k-point mesh directly with respect to the coordinate axes in reciprocal space, rather then with respect to the reciprocal lattice vectors. This can be achieved by setting t_kpoint_rel to .false.. The unit of length on the coordinate axes is $2\pi/a_{\rm lat}$ in this case. The folding parameters can be used as well to enhance the number of k-points, if desired. However, one should keep in mind that the k-point sets specified in that way might have little symmetry, i.e. their number is not significantly reduced by the built-in symmetry reduction algorithm of fhi98start.

• Reduced k-points and symmetry
  Apart from the translational symmetry of the Bravais lattice, the crystal structure under investigation may often have additional point group symmetries. These can be used to reduce the number of k-points which are needed in the actual calculation (and thus the memory demand) substantially. To perform the integrals in the Brillouin zone , it is sufficient to sample the contribution from a subset of non-symmetry-equivalent k-points only. Therefore the integrand (e.g. the charge density) is calculated only at these points. The integrand with the full symmetry can be recovered from its representation by non-symmetry-equivalent k-points whenever this is required.  

  The fhi98start utility is set up to automatically exploit these point group symmetries. First, the point group symmetry operations applicable to the unit cell are determined and stored in the form of symmetry matrices. Secondly, fhi98start seeks to reduce the elements of the k-point mesh to the q u subset which is irreducible under those symmetry matrices. Only this subset is forwarded in the inp.ini  file for further use in the main computations. The performance of the reduction procedure can be monitored by inspecting the output in the file start.out . An estimate for the sampling quality of the k-point set is given on the basis of the analysis of 'shells' (see Chadi and Cohen[31] for details). For a good k-point set, the contribution from the leading 'shells' should vanish. Some comments for interested users:

  • For a slab k-point set, those shells that contain contributions from lattice vectors with a finite z-component cannot vanish, thus they must be disregarded when judging the quality of the basis set.
  • The quality assessment only makes sense for systems with a band gap. The effect of having a sharp integration boundary, the Fermi surface, for a metal is not accounted for by Chadi and Cohen's shell analysis.

  Note: Even if there are no point group symmetries, the vectors k and $-{\bf k}$ are symmetry-equivalent by virtue of time-inversion symmetry. For this reason, the number of k-points is reduced by at least a factor of two for any reasonable choice of a k-point mesh.