Setting up k-points and spin in ONETEP

Author:: Chengcheng Xiao, Miguel Escobar Azor and Nicholas D.M. Hine
Date:: May 2025

Spin-polarisation in ONETEP

TBC

k-points and Brillouin zone sampling in ONETEP

ONETEP now supports k-point sampling in the Brillouin zone. By utilising the k-point sampling technique, the periodicity of the system can be exploited to reduce the computational cost of the calculation and small systems (from unit cells upwards) can be efficiently simulated.

The k-point sampling technique is implemented in two modes: the plane-wave (PW) mode and the tight-binding (TB) mode. The default behaviour in ONETEP is to assume the system is periodic but that the simulation cell is large enough that k-point sampling is not required. In order to launch a calculation on a simulation cell where this is not the case, so k-point sampling is required, the user must first choose between these modes as appropriate to the periodicity of their system, and then specify a grid of k-points on which the calculation is to be performed.

Plane-wave mode (PW) mode

In the plane-wave (PW) mode, the k-dependent Bloch functions are expressed using cell-periodic functions and and the Bloch phase factor. The cell-periodic functions are further expanded in terms of NGWFs.

Specifically, to satisfy the minimum image convention (i.e., no NGWF should overlap with its periodic images), the localisation constraint is relaxed along the directions where k-point sampling is applied. Since the NGWFs are expressed using psinc functions, by relaxing the localisation constraint, NGWFs are essentially expressed using plane waves (hence the name “plane-wave mode”).

Theory

The Schrödinger equation (or in the case of DFT, the Kohn-Sham equation) for a Bloch eigenstate \(\psi_{n,\mathbf{k}}\) is given by

(2)\[\hat{H} \psi_{n,\mathbf k} (\mathbf r) = \epsilon_{n,\mathbf k} \psi_{n,\mathbf k} (\mathbf r).\]

where \(\mathbf k\) labels the k-point, \(n\) the band number. Immediately we can see that we need to solve this equation separately at each k-point, hence we need to have seperate representations of any quantities arising from \(\psi_{n,\mathbf k}\) for each k-point.

The most generic form of the Bloch functions is written as,

(3)\[\psi_{n,\mathbf k} (\mathbf r) = e^{i\mathbf k \cdot \mathbf r} u_{n,\mathbf k} (\mathbf r),\]

where \(u_{n,\mathbf k} (\mathbf r)\) is a function with the periodicity of the unit cell, i.e.,

\[u_{n,\mathbf k} (\mathbf r) = u_{n,\mathbf k} (\mathbf r + \mathbf R),\]

and \(\mathbf R\) is a Bravais lattice vector.

Inserting (3) into (2), we get

\[\hat{H} e^{i\mathbf k \cdot \mathbf r} u_{n,\mathbf k} (\mathbf r) = \epsilon_{n,\mathbf k} e^{i\mathbf k \cdot \mathbf r} u_{n,\mathbf k} (\mathbf r) .\]

Premultiplying both sides by \(e^{i\mathbf k \cdot \mathbf{r}}\), we get

(4)\[e^{-i\mathbf k \cdot \mathbf r} \hat{H} e^{i\mathbf k \cdot \mathbf r} u_{n,\mathbf k} (\mathbf r) = \epsilon_{n,\mathbf k} u_{n,\mathbf k} (\mathbf r).\]

If we introduce the k-dependent Hamiltonian operator \(\hat{H}(\mathbf k)\),

\[\hat{H}(\mathbf k) = e^{-i\mathbf k \cdot \mathbf r} \hat{H} e^{i\mathbf k \cdot \mathbf r},\]

(4) can then be written as

(5)\[\hat{H}(\mathbf k) u_{n,\mathbf k} (\mathbf r) = \epsilon_{n,\mathbf k} u_{n,\mathbf k} (\mathbf r).\]

In ONETEP we use a set of k-dependent NGWFs (\(\{ \phi_\alpha^{\mathbf k}\}\)) as our support functions, so the cell-periodic part of the Bloch wavefunction \(u_{n,\mathbf{k}}\) can be written as

\[u_{n,\mathbf k} (\mathbf r) = \sum_{\alpha} c_{n,\mathbf k}^{\alpha} \phi_{\alpha}^{\mathbf k} (\mathbf r),\]

where \(c_{n,\mathbf k}^{\alpha}\) is a non-unitary rotation matrix that rotates the NGWFs into eigenstates.

Expanding \(u_{n,\mathbf k}\) using in terms of NGWFs (which, as we will see later will need to be “extended” along periodic directions in PW mode), (5) can be expressed as

(6)\[\hat{H}(\mathbf k) u_{n,\mathbf k} (\mathbf r)= \sum_{\beta} c_{n, \mathbf k}^{\beta}\hat{H}(\mathbf k) \phi_{\beta}^{\mathbf k} (\mathbf r) = \epsilon_{n,\mathbf k} \sum_{\beta} c_{n,\mathbf k}^{\alpha} \phi_{\beta}^{\mathbf k} (\mathbf r)\]

If we multiply both sides of (6) by \(\phi_{\alpha}^{\mathbf k *}(\mathbf r)\), and integrate over \(\mathbf r\), we get

(7)\[\sum_{\beta} c_{n, \mathbf k}^{\beta} \int \phi_{\alpha}^{\mathbf k *}(\mathbf r)\hat{H}(\mathbf k) \phi_{\beta}^{\mathbf k} (\mathbf r) d\mathbf r= \sum_{\beta} c_{n, \mathbf k}^{\beta} \epsilon_{n,\mathbf k} \int \phi_{\alpha}^{\mathbf k *}(\mathbf r) \phi_{\beta}^{\mathbf k} (\mathbf r) d\mathbf r.\]

In matrix notation, (7) can be written as

(8)\[\sum_{\beta} H_{\alpha \beta}(\mathbf k) c_{n, \mathbf k}^{\beta} = \epsilon_{n,\mathbf k} \sum_{\beta} c_{n, \mathbf k}^{\beta} S_{\alpha \beta}^{\mathbf k},\]

where the Hamiltonian matrix \(H\) elements are:

\[H_{\alpha \beta} (\mathbf k) = \int \phi_{\alpha}^{\mathbf k *}(\mathbf r)\hat{H}(\mathbf k) \phi_{\beta}^{\mathbf k} (\mathbf r) d\mathbf r.\]

and the overlap matrix \(S\) elements are:

\[S_{\alpha \beta} (\mathbf k) = \int \phi_{\alpha}^{\mathbf k *}(\mathbf r)\phi_{\beta}^{\mathbf k} (\mathbf r) d\mathbf r.\]

(8) is the generalized Kohn-Sham equation under the basis of \(\phi\). If \(\phi\) were simply plane waves, \(\mathbf{S}\) would be the identity, and we could return to the canonical expression of the Kohn-Sham equation in the plane-wave basis:

\[\sum_{\beta} H_{\alpha \beta}(\mathbf k) c_{n, \mathbf k}^{\beta} = \epsilon_{n,\mathbf k} c_{n, \mathbf k}^{\alpha}\]

It is worth nothing that in the plane-wave basis the explicit k-dependence only exists in the Hamiltonian matrix and the k-index of eigenvectors is the result of solving such a k-dependent Hamiltonian. In the NGWF basis, there is also k-dependence in the overlap matrix since it is no longer identity.

The Hamiltonian can be expanded into three terms: the kinetic energy term [\(T_{\alpha \beta}(\mathbf{k})\)], the potential term (including local potential from atomic cores and exchange correlation terms) and the non-local terms from pseudopotentials and we will now see how these terms are expressed with k-dependence.

Kinetic energy term

The matrix elements of the kinetic energy operator under the psinc basis are given by

\[T_{\alpha \beta}(\mathbf{k}) =\int d^3 r \phi_\alpha^{\mathbf{k} *}(\mathbf{r}) \mathrm{e}^{-\mathrm{i} \mathbf{k} \cdot \mathbf{r}} (-\frac{\nabla_{\mathbf r}^2}{2}) \mathrm{e}^{\mathrm{i} \mathbf{k} \cdot \mathbf{r}} \phi_\beta^{\mathbf{k}}(\mathbf{r}),\]

Using psinc functions \(D_i(\mathbf r)\):

\[D_{i} (\mathbf r) = \frac{1}{N}\sum_p e^{i\mathbf G_p \cdot (\mathbf r - \mathbf r_i)},\]

The kinetic matrix elements becomes

\[\begin{split}\begin{aligned} T_{\alpha \beta}(\mathbf{k}) &=\int d^3 r \phi_\alpha^{\mathbf{k} *}(\mathbf{r}) \mathrm{e}^{-\mathrm{i} \mathbf{k} \cdot \mathbf{r}} (-\frac{\nabla_{\mathbf r}^2}{2}) \mathrm{e}^{\mathrm{i} \mathbf{k} \cdot \mathbf{r}} \phi_\beta^{\mathbf{k}}(\mathbf{r}) \\ &= -\frac{V}{2N^2} \sum_{ij} c^{\mathbf k *}_{j \alpha} c^{\mathbf k}_{i \beta} \sum_{p} (-|\mathbf G_p + \mathbf k|^2) e^{i\mathbf G_p (\mathbf r_j - \mathbf r_i)}. \\ &= \frac{V}{2N^2} \sum_{ij} c^{\mathbf k *}_{j \alpha} c^{\mathbf k}_{i \beta} \sum_{p} (|\mathbf G_p + \mathbf k|^2) e^{i\mathbf G_p (\mathbf r_j - \mathbf r_i)}. \end{aligned}\end{split}\]

where we see that we only need to add k vector to the \(\mathbf{G}\) vectors to make the kinetic energy matrix k-dependent.

Local potential, hartree and exchange-correlation energy terms

The k-dependent local potential energy matrix is given by:

\[\begin{split}\begin{aligned} V_{\mathrm{LHXC},\alpha\beta}(\mathbf{k}) &= \int d^3 r \phi_\alpha^{\mathbf{k} *}(\mathbf{r}) e^{-\mathrm{i} \mathbf{k} \cdot \mathbf{r}} V_\mathrm{LHXC}(\mathbf{r}) e^{\mathrm{i} \mathbf{k} \cdot \mathbf{r}} \phi_\beta^{\mathbf{k}}(\mathbf{r}) \\ &= \int d^3 r \phi_\alpha^{\mathbf{k} *}(\mathbf{r}) V_\mathrm{LHXC}(\mathbf{r}) \phi_\beta^{\mathbf{k}}(\mathbf{r}) \\ \end{aligned}\end{split}\]

And we see that no explicit k-dependence is needed in the local potential term.

Non-local potential term

The k-dependent non-local potential term is given by:

\[V_{\mathrm{NL},\alpha\beta}(\mathbf{k}) = \sum_{I,lm} \left<\phi_\alpha^{\mathbf{k}}|e^{-\mathrm{i} \mathbf{k} \cdot \mathbf{r}}|\xi^I_{lm}\right> V_{lm}^I(\mathbf{k}) \left<\xi^I_{lm}|e^{\mathrm{i} \mathbf{k} \cdot \mathbf{r}}|\phi_\beta^{\mathbf{k}}\right>\]

and we see that the phase factors are to be augmented to the non-local projectors so that the non-local potential term is k-dependent.

Tight-binding (TB) mode

The tight-binding (TB) mode is designed to use fully localised NGWFs ( extended_ngwf : F F F) and the k-point sampling is performed by augmenting the Hamiltonian, overlap and other matrices with k-dependent phase factors.

Specifically, in the TB mode, we adopt the Bloch sum form of the Bloch function:

\[\psi_{n\mathbf k}(\mathbf r) = \sum_{\mathbf R} e^{i\mathbf k \cdot \mathbf R} \sum_\alpha c_{n,\mathbf k}^\alpha \phi_\alpha(\mathbf r - \mathbf R)\]

where \(c_{n\mathbf k}\) is the k-dependent coefficient that rotates the basis into Kohn-Sham eigenstates (similar to that in the PW mode) and \(\mathbf{R}\) is the lattice vector of the unit cell. Using this expression, we can re-express the charge density and the Kohn-Sham energy with k-dependent phase factors.

Using this Bloch sum, the charge density becomes

\[\begin{split}\begin{aligned} n(\mathbf r) =& \sum_n f_{n,\mathbf k} \sum_{\mathbf k} w_{\mathbf k} |\psi_{n,\mathbf k}(\mathbf r)|^2\\ =&N_\mathrm{cell}\sum_{\alpha\beta}\sum_{\mathbf k} w_{\mathbf k} K_{\mathbf k}^{\beta\alpha} \sum_{\mathbf R'} e^{i\mathbf k\cdot \mathbf R'} \phi_{\alpha}(\mathbf r-\mathbf 0) \phi_{\beta}(\mathbf r - \mathbf R') \end{aligned}\end{split}\]

where \(f_{n\mathbf k}\) is the occupation of the Bloch state at band number \(n\) and \(\mathbf k\), \(w_{\mathbf k}\) is the k-point weight. The density kernel element \(K^{\alpha\beta}_ {\mathbf k}\) is expressed as \(\sum_{n} c_{n\mathbf k}^{\alpha*} f_{n\mathbf k} c_{n\mathbf k}^\beta\) and is k-dependent.

Since NGWFs are localised and only overlap with certain other NGWFs, we can limit the summation over \(\mathbf R'\) over all cells to only the one that makes NGWF \(\alpha\) overlaps with \(\beta\). Doing this with the phase factors written in a product form gives

\[\begin{aligned} n(\mathbf r) &=N_\mathrm{cell}\sum_{\alpha\beta}\sum_{\mathbf k} w_{\mathbf k} K_{\mathbf k}^{\beta\alpha} \phi_{\alpha}(\mathbf r - \mathbf 0) \phi_{\beta}(\mathbf r - \mathbf R') \prod_i \theta(k_i,r_{\alpha i} - r_{\beta i}, R'_i) \end{aligned}\]

where \(\mathbf R'\) moves the \(\beta\)-NGWF to be the nearest-neighbor of the \(\alpha\)-NGWF based on the relative location of their centers \(r_{\alpha}\) and \(r_{\beta}\). The one-dimensional phase factor \(\theta\) takes the form along the \(i\)-th lattice coordinate as

\[\begin{split}\theta(k, r, R)= \begin{cases}1 & |r| \leqslant \frac{1}{2} \\ e^{i k R} & r>\frac{1}{2} \\ e^{-i k R} & r<-\frac{1}{2}\end{cases}.\end{split}\]

Or, in a shorthand notation we can write it as \(\Theta[\mathbf{k},\alpha,\beta]\).

Similarly, the Hamiltonian matrix elelments can be expressed as

\[H_{\mathbf{k},\alpha\beta} = H_{\alpha\beta} \Theta[\mathbf{k},\alpha,\beta]\]

Special treatment is needed for the non-local pseudopotential terms where phase factors are added based on the relative location between NGWFs and the non-local projectors. Moreover, the gradient contribution also needs to be augmented with phase factors.

One thing to note is that when performing tensor corrections, the k-independent overlap matrix is used instead of the one augmented with phase factors, this is due to the fact that only one set of NGWFs is used for all k-points in the TB mode.

Brillouin zone sampling

To find the ground state of the system, we need to sample the Brillouin zone by performing a summation over the results of different k-points. Specifically, we need to perform k-point average over all energy components and the charge density.

The k-point sampling is performed using the Monkhorst-Pack scheme where kpoints along a spcific direction are generated (in the first Brillouin zone) by

\[\frac{2r-q-1}{2q}\]

where \(q\) is the number of k-points in each direction and \(r\) is the k-point index. An optional shift of half a grid cell can be added so that \(\Gamma\) point is included in the sampling.

There are two ways to define the k-point sampling in ONETEP:

Automatic generation using the Monkhorst-Pack scheme.:
```
kpoint_grid_shift : 0 1 1
kpoint_grid_size : 3 6 6
```
which indicates that the size of the k-point grid is 3 along a direction, 6 along b direction (shifted by half a grid distance) and 6 along c direction (also shifted by half a grid distance).

Manual generation using the k-point coordinates:

%block kpoints_list
    0.5000000000    0.50000000    0.0000000000    0.125
    0.5000000000   -0.50000000    0.0000000000    0.125
   -0.5000000000    0.50000000    0.0000000000    0.125
   -0.5000000000   -0.50000000    0.0000000000    0.125
    0.5000000000    0.50000000    0.5000000000    0.125
    0.5000000000   -0.50000000    0.5000000000    0.125
   -0.5000000000    0.50000000    0.5000000000    0.125
   -0.5000000000   -0.50000000    0.5000000000    0.125
%endblock kpoints_list

where k-points are defined in the fractional coordinates and the last column is the weight of the k-point.

Kpar parallelisation

Since the calculations with k-points are almost fully isolated with one and another, it makes sense to use k-pool parallelisation where multiple instances of ONETEP (k-parallelisation groups, or kpars in short) is launched altogether, each in charge of running a sub-set of k-points (and in each kpar, k-points are looped over in serial).

The kpar parallelisation is controlled by the keyword num_kpar in the input file:

num_kpars : 4

This means that the k-points are divided into 4 groups and 4 ONETEP instances will be launched.

When kpar parallelisation is used, the total number of MPI processes will be divided by the number of kpars, i.e., if you have 16 MPI processes and 4 kpars, each kpar will use 4 MPI processes. This means that the number of processes needs to be divisible by the number of kpars. Also, because each ONETEP instance can only use the number of processes less or equal to the number of atoms in the system, the total number of MPI processes should not exceed number of kpars times the number of atoms in the system.

OpenMP threading can also be used in combination with MPI parallelisation and kpar parallelisation. Since one standalone instance of ONETEP is launched for each kpar, the number of OpenMP threads only takes effect within each kpar, attaching to each ONETEP instance. This means that keywords such as THREADS_NUM_FFTBOXES and THREADS_MAX (BASIC) will only affect threads within each instance of ONETEP.

To give an example, one can run a calculation with: 4 kpars, 4 MPI processes per kpar and 4 OpenMP threads per MPI process.

This means one would require a total number of 64 cores (1 node), used by 16 MPI processes and 4 openMP threads per MPI process. An example PBS script for this job is given below:

#PBS -l select=1:ncpus=64:mpiprocs=16:ompthreads=4:mem=200gb

Hybrid and extended NGWFs

In the PW mode, NGWFs needs to be extended. This is turned on by using the keyword extended_ngwf in the input file. e.g.,

extended_ngwf : T T T

It is also possible to only allow NGWFs to be extended along certain directions, hence utlising the periodicity of the system only along those directions via k-point sampling.

It is worth nothing that the extended NGWFs can also run withotu k-point sampling and should produce the same results as running a plan-wave DFT calculation with only one k-point (i.e., the \(\Gamma\) point).

Fixed kernel calculation

When using fully extended NGWFs (extended_ngwf : T T T) it is possible to perform fixed kernel calculation where the number of NGWFs used is the same as the number of occupied states. Since we are not optimising the density kernel, this will only work with known insulators (i.e., all NGWFs are fully occupied).

This is done by setting:

maxit_lnv=-1
minit_lnv=-1

One thing to note is that to use this feature, sometimes one has to turn on permit_unusual_ngwf_count in the input file and play around with the number of NGWFs for each species so that the total number of NGWFs is equal to the number of occupied states.

Additional notes

Currently, full Brillouin zone sampling is only tested for norm-conserving pseudopotentials. Here’s a brief list of supported functionalities:

Ground state energy calculation with LNV (exact_lnv : T) and EDFT (edft : T).
Geometry optimisation (but no cell-optimisation).
Parts of the properties module (e.g., charge density outputs, eigenvalue outputs).

Keywords

extended_ngwf [Basic, bool bool bool, default F F F]. Turn on extended NGWFs along the three directions.
kpoint_method [Basic, default None]. The method used to generate the k-point grid. The options are:
- PW: Plane-wave mode. Requires NGWFs to be extended along the periodic directions where k-point sampling is applied.
- TB: Tight-Binding mode. Requires NGWFs to be fully localised.
- None: No k-point sampling.
kpoint_grid_shift [Basic int int int, default 0 0 0]. The shift of the k-point grid.
kpoint_grid_size [Basic int int int, default 1 1 1]. The size of the k-point grid.
num_kpars [Basic int, default 1]. The number of k-parallelisation groups.