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 points. The details of this procedure are as follows: In**fhi98md**, the Brillouin zone is spanned by the reciprocal lattice vectors and 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 and 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.

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 likeParameter 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 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 -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 , the other with a unit cell. In this case, useParameter Value *i_facs(1..3)*`4 8 1`**k**-point folding factorsin the first case, and

Parameter Value *i_facs(1..3)*`2 4 1`**k**-point folding factorsin the second case, leaving the other parameters unchanged.

**Note:**is orthogonal to the real lattice vectors and . If is the long edge of your real space unit cell, spans the short edge of your Brillouin zone. Therefore, the**k**-point sampling mesh has fewer points in the direction and more points in the 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 followsParameter 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 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 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.- For a slab