Tutorial kpoint convergence testing
Author: Jane E. Herriman
You may want to take a look at the ground state [tutorial] (https://github.com/LLNL/qball/wiki/Tutorial-Ground-State) to familiarize yourself with Qb@ll input variables before doing this tutorial.
In this tutorial, we will see how to do kpoint convergence testing. The system with which we will be working is cubic, so our mesh will have the same number of kpoints along each reciprocal space lattice vector. That means we want to find the smallest value N such that the calculated energy of our system is not significantly different for a (N x N x N) mesh and a (M x M x M) mesh where M > N.
We've pre-generated some equipartitioned kpoint meshes for you to work with. You can find them in the Appendix below, along with an example of how to pick a kpoint mesh for a non-cubic system.
Create an empty text file to use as your input file. We suggest you call it GaAs.i. Add the variables/commands discussed in the next two sections to this input file.
We will be using a primitive cell of GaAs in this tutorial. The commands below provide information about this system, such as the geometry of a GaAs primitive cell, and should be included in your input file.
You will want to define the lattice vectors of your cell. Recall that units are in Bohr radii (a.u.).
set cell 0.00000 3.30000 3.30000 3.30000 0.00000 3.30000 3.30000 3.30000 0.0000
Next, define the atomic species present in your cell -- create an association between an xml file containing all information that describes this particular species (including pseudopotential information) and an alias for the species. The xml files listed below can be found [here] (https://raw.githubusercontent.com/LLNL/qball/master/examples/tutorial/pseudos/As_HSCV_PBE-1.0.xml) and here.
species Gallium Ga_HSCV_PBE-1.0.xml
species Arsenic As_HSCV_PBE-1.0.xml
Define the atoms in the unit cell and their positions (in a.u.).
atom Ga1 Gallium 0.0000 0.0000 0.0000
atom As1 Arsenic 1.6500 1.6500 1.6500
Make sure to set variables discussed in the ground state tutorial, namely ecut, ecutprec, wf_dyn, and threshold_scf. Try using values for ecut and ecutprec of 45.0 and 8.0 Rydberg, respectively.
To do kpoint convergence testing, you'll have to run several calculations, each using a different mesh on which kpoints are defined. For each run, you'll want to set the total number of kpoints in your input file via set nkpoints. Note that this is the number of irreducible kpoints; for example, a 2x2x2 mesh has 8 points, but only 4 irreducible kpoints by symmetry. Next, you will have to define the coordinates and weight for each kpoint with the kpoint command. You can copy the commands necessary to set up a (N x N x N) kpoint mesh for N between 1 and 7 into your input file from the Appendix below.
Finally, you'll want to add the randomize_wf and run commands to the end of your input file. You may want to set the maximum number of self-consistent iterations to a few hundred, e.g. run 0 400.
After setting up your input file and adding the requisite *.xml files to your working directory, you should be ready to run qball. You want to catch the output in some output file, e.g. GaAs.out via
qball GaAs.i > GaAs.out
assuming qball is the path to your qball binary and that your input file is named GaAs.i. Once this calculation is complete, search for the last value of etotal printed to your output file via
grep "<etotal>" GaAs.out | tail -1
Make a note of this value for the total energy and the kpoint mesh that you used to generate this value. Run similar versions of this calculation with different kpoint meshes and record the resulting values of the total energies.
When you have run several calculations with different kpoint meshes, plot the energy of the GaAs unit cell as a function of the number of kpoints used. Using kpoint (NxNxN) kpoint meshes for N between 1 and 10, we got the following plot: 
but we can see a bit more if we zoom in: 
.
You can see from the images linked above that for equipartitioned kpoint meshes with more than perhaps 27 kpoints (a 3x3x3 mesh), the energy of the GaAs unit cell seems to oscillate. The amplitude of the oscillations should decrease with the number of kpoints used. The number of kpoints that you use should therefore be determined by your application and the accuracy required.
Below are provided lists of coordinates that define equipartitioned meshes in reciprocal space for a primitive cell of zinc blende GaAs. Each coordinate comes with an associated weight.
Coordinates for a 1x1x1 mesh:
set nkpoints 1
kpoint 0.00000000000 0.00000000000 0.00000000000 1.0000000000000
Coordinates for a 2x2x2 mesh:
set nkpoints 4
kpoint -0.250000 -0.250000 -0.250000 0.125000
kpoint -0.250000 -0.250000 0.250000 0.375000
kpoint -0.250000 0.250000 0.250000 0.375000
kpoint 0.250000 0.250000 0.250000 0.125000
Coordinates for a 3x3x3 mesh:
set nkpoints 5
kpoint -0.333333 -0.333333 -0.333333 0.148148
kpoint -0.333333 -0.333333 0.000000 0.222222
kpoint -0.333333 -0.333333 0.333333 0.444444
kpoint -0.333333 0.000000 0.000000 0.148148
kpoint 0.000000 0.000000 0.000000 0.037037
Coordinates for a 4x4x4 mesh:
set nkpoints 20
kpoint -0.375000 -0.375000 -0.375000 0.015625
kpoint -0.375000 -0.375000 -0.125000 0.046875
kpoint -0.375000 -0.375000 0.125000 0.046875
kpoint -0.375000 -0.375000 0.375000 0.046875
kpoint -0.375000 -0.125000 -0.125000 0.046875
kpoint -0.375000 -0.125000 0.125000 0.093750
kpoint -0.375000 -0.125000 0.375000 0.093750
kpoint -0.375000 0.125000 0.125000 0.046875
kpoint -0.375000 0.125000 0.375000 0.093750
kpoint -0.375000 0.375000 0.375000 0.046875
kpoint -0.125000 -0.125000 -0.125000 0.015625
kpoint -0.125000 -0.125000 0.125000 0.046875
kpoint -0.125000 -0.125000 0.375000 0.046875
kpoint -0.125000 0.125000 0.125000 0.046875
kpoint -0.125000 0.125000 0.375000 0.093750
kpoint -0.125000 0.375000 0.375000 0.046875
kpoint 0.125000 0.125000 0.125000 0.015625
kpoint 0.125000 0.125000 0.375000 0.046875
kpoint 0.125000 0.375000 0.375000 0.046875
kpoint 0.375000 0.375000 0.375000 0.015625
Coordinates for a 5x5x5 mesh:
set nkpoints 14
kpoint -0.400000 -0.400000 -0.400000 0.032000
kpoint -0.400000 -0.400000 -0.200000 0.096000
kpoint -0.400000 -0.400000 0.000000 0.048000
kpoint -0.400000 -0.400000 0.200000 0.096000
kpoint -0.400000 -0.400000 0.400000 0.096000
kpoint -0.400000 -0.200000 -0.200000 0.096000
kpoint -0.400000 -0.200000 0.000000 0.096000
kpoint -0.400000 -0.200000 0.200000 0.192000
kpoint -0.400000 0.000000 0.000000 0.032000
kpoint -0.400000 0.000000 0.200000 0.096000
kpoint -0.200000 -0.200000 -0.200000 0.032000
kpoint -0.200000 -0.200000 0.000000 0.048000
kpoint -0.200000 0.000000 0.000000 0.032000
kpoint 0.000000 0.000000 0.000000 0.008000
Coordinates for a 6x6x6 mesh:
set nkpoints 56
kpoint -0.416667 -0.416667 -0.416667 0.004630
kpoint -0.416667 -0.416667 -0.250000 0.013889
kpoint -0.416667 -0.416667 -0.083333 0.013889
kpoint -0.416667 -0.416667 0.083333 0.013889
kpoint -0.416667 -0.416667 0.250000 0.013889
kpoint -0.416667 -0.416667 0.416667 0.013889
kpoint -0.416667 -0.250000 -0.250000 0.013889
kpoint -0.416667 -0.250000 -0.083333 0.027778
kpoint -0.416667 -0.250000 0.083333 0.027778
kpoint -0.416667 -0.250000 0.250000 0.027778
kpoint -0.416667 -0.250000 0.416667 0.027778
kpoint -0.416667 -0.083333 -0.083333 0.013889
kpoint -0.416667 -0.083333 0.083333 0.027778
kpoint -0.416667 -0.083333 0.250000 0.027778
kpoint -0.416667 -0.083333 0.416667 0.027778
kpoint -0.416667 0.083333 0.083333 0.013889
kpoint -0.416667 0.083333 0.250000 0.027778
kpoint -0.416667 0.083333 0.416667 0.027778
kpoint -0.416667 0.250000 0.250000 0.013889
kpoint -0.416667 0.250000 0.416667 0.027778
kpoint -0.416667 0.416667 0.416667 0.013889
kpoint -0.250000 -0.250000 -0.250000 0.004630
kpoint -0.250000 -0.250000 -0.083333 0.013889
kpoint -0.250000 -0.250000 0.083333 0.013889
kpoint -0.250000 -0.250000 0.250000 0.013889
kpoint -0.250000 -0.250000 0.416667 0.013889
kpoint -0.250000 -0.083333 -0.083333 0.013889
kpoint -0.250000 -0.083333 0.083333 0.027778
kpoint -0.250000 -0.083333 0.250000 0.027778
kpoint -0.250000 -0.083333 0.416667 0.027778
kpoint -0.250000 0.083333 0.083333 0.013889
kpoint -0.250000 0.083333 0.250000 0.027778
kpoint -0.250000 0.083333 0.416667 0.027778
kpoint -0.250000 0.250000 0.250000 0.013889
kpoint -0.250000 0.250000 0.416667 0.027778
kpoint -0.250000 0.416667 0.416667 0.013889
kpoint -0.083333 -0.083333 -0.083333 0.004630
kpoint -0.083333 -0.083333 0.083333 0.013889
kpoint -0.083333 -0.083333 0.250000 0.013889
kpoint -0.083333 -0.083333 0.416667 0.013889
kpoint -0.083333 0.083333 0.083333 0.013889
kpoint -0.083333 0.083333 0.250000 0.027778
kpoint -0.083333 0.083333 0.416667 0.027778
kpoint -0.083333 0.250000 0.250000 0.013889
kpoint -0.083333 0.250000 0.416667 0.027778
kpoint -0.083333 0.416667 0.416667 0.013889
kpoint 0.083333 0.083333 0.083333 0.004630
kpoint 0.083333 0.083333 0.250000 0.013889
kpoint 0.083333 0.083333 0.416667 0.013889
kpoint 0.083333 0.250000 0.250000 0.013889
kpoint 0.083333 0.250000 0.416667 0.027778
kpoint 0.083333 0.416667 0.416667 0.013889
kpoint 0.250000 0.250000 0.250000 0.004630
kpoint 0.250000 0.250000 0.416667 0.013889
kpoint 0.250000 0.416667 0.416667 0.013889
kpoint 0.416667 0.416667 0.416667 0.004630
Coordinates for a 7x7x7 mesh:
set nkpoints 30
kpoint -0.428571 -0.428571 -0.428571 0.011662
kpoint -0.428571 -0.428571 -0.285714 0.034985
kpoint -0.428571 -0.428571 -0.142857 0.034985
kpoint -0.428571 -0.428571 0.000000 0.017493
kpoint -0.428571 -0.428571 0.142857 0.034985
kpoint -0.428571 -0.428571 0.285714 0.034985
kpoint -0.428571 -0.428571 0.428571 0.034985
kpoint -0.428571 -0.285714 -0.285714 0.034985
kpoint -0.428571 -0.285714 -0.142857 0.069971
kpoint -0.428571 -0.285714 0.000000 0.034985
kpoint -0.428571 -0.285714 0.142857 0.069971
kpoint -0.428571 -0.285714 0.285714 0.069971
kpoint -0.428571 -0.285714 0.428571 0.069971
kpoint -0.428571 -0.142857 -0.142857 0.034985
kpoint -0.428571 -0.142857 0.000000 0.034985
kpoint -0.428571 -0.142857 0.142857 0.069971
kpoint -0.428571 0.000000 0.000000 0.011662
kpoint -0.428571 0.000000 0.142857 0.034985
kpoint -0.428571 0.000000 0.285714 0.034985
kpoint -0.428571 0.285714 0.285714 0.034985
kpoint -0.285714 -0.285714 -0.285714 0.011662
kpoint -0.285714 -0.285714 -0.142857 0.034985
kpoint -0.285714 -0.285714 0.000000 0.017493
kpoint -0.285714 -0.142857 -0.142857 0.034985
kpoint -0.285714 -0.142857 0.000000 0.034985
kpoint -0.285714 0.000000 0.000000 0.011662
kpoint -0.142857 -0.142857 -0.142857 0.011662
kpoint -0.142857 -0.142857 0.000000 0.017493
kpoint -0.142857 0.000000 0.000000 0.011662
kpoint 0.000000 0.000000 0.000000 0.002915
The densities of kpoints along each reciprocal space basis vector should be equivalent, or approximately equivalent. When working with a cubic lattice, this is relatively straight-forward: we can simply use an equipartitioned (N x N x N) mesh in reciprocal space.
How should our mesh of kpoints change for a lower symmetry lattice, i.e. one that is non-cubic? Let's call the basis vectors used to describe our lattice are a, b, and c. In reciprocal space, the lengths of the corresponding vectors are |1⁄a|, |1⁄b|, and |1⁄c|. To ensure equal densities of kpoints along each reciprocal space basis vector, we want to choose a (U x V x W) mesh in reciprocal space such that |U⁄a| == |V⁄b| == |W⁄c|.
Here's a quick example: Let's say we're working with a unit cell of an orthorhombic lattice. In this particular orthorhombic lattice, b == c == 2a. Substituting into one of the equations above, |U⁄a| == |V⁄2a| == |W⁄2a|. This means that 2U == V == W. You'll therefore want to use reciprocal space meshes of the form (U x 2U x 2U) for this system.
This should by no means be taken as an introduction to reciprocal/fourier space. If you're unfamiliar with the above concepts, you may want to check out an introductory text in solid state physics such as Kittel.