Danny Vanpoucke

One of the main differences between theory and computational research is the fact that the latter has to deal with finite resources; mainly time and storage. Where theoretical calculations involve integrations over continuous spaces, infinite sums and basis sets, computational work performs numerical integrations as weighted sums over finite grids and cuts of infinite series. As an infinite number of operations would take an infinite amount of time, it is clear why numerical evaluations are truncated. If the contributions of an infinite series become smaller and smaller, it is also clear that at some point the contributions will become smaller than the numerical accuracy, so continuation beyond that point is …pointless.

In case of ab initio quantum mechanical calculations, we aim to get as accurate results at an as low computational cost as possible. Even with the current availability of computational resources, an infinite sum would still take an infinite amount of time. In addition, although parallelization can help out a lot in getting access to additional computational resources during the same amount of real time, codes are not infinitely parallel, so at some point adding more CPU’s will no longer speed up the calculations. Two important parameters in quantum mechanical calculations to play with are the basis set size (Or kinetic energy cut off, in case of plane wave basis sets. In which case this can also be related to the real space integration grid) and the integration grid for the reciprocal space (the k-point grid).

These two parameters are not unique to VASP, they are present in all quantum mechanical codes, but we will use VASP as an example here. The example system we will use is the α-phase of Cerium, using the PBE functional. The default cut-off energy used by VASP is 299 eV.

1.     Basis set size/Kinetic energy cut-off

What a basis set is and how it is defined depends strongly on the code. As such you are referred to the manual/tutorials of your code of interest.(VASP workshop) One important thing to remember, however, is the fact that although a plane wave basis set is “nicely behaved” (bigger basis = more accurate result) this is not true for all types of basis sets (Gaussian basis sets are an important example here).

How do you perform a convergence test?

1. Get a geometry of your system of interest.

   This does not need to be a fully optimized geometry, an experimental geometry or a reasonable manually constructed geometry will do fine, as long as it gives you a converged result at the end of your static calculation. A convergence test should not depend on the exact geometry of your system. Rather it should tell you how well your setting converges your result with regard to the energy found on the potential energy surface.

2. Fix all other settings

   (to reasonable values, although the settings should—to make your life somewhat sane—be independent with regard to convergence testing).
   VASP specific parameters of importance:

   • PREC : should be at least normal, but high or accurate are also possible
   • EDIFF : a value of 1.0E-6 to 1.0E-8 are reasonable for small systems. Note that this value should be much smaller than the accuracy you wish to obtain.
   • NSW = 0; IBRION = -1 : It should be static calculations.
   • ISPIN : If you intend to perform spin polarized calculations, you should also include this in your convergence tests. Yes it increases the computational cost, but remember that convergence tests will only take a fraction of the computational costs of your project, and can save you a lot of work and resources later on.
   • NBANDS : You may want to manually fix the number of electronic bands, which will allow for comparison of timing results.
   • LCHARG = .TRUE. ; ICHARG = 1: If you are not that interested in timing (or use average time of electronic loops instead of total CPU time), and want to speed things up a bit, you can use the electron density from a cheaper calculation as a starting point.
   • KPOINTS-file: use a non-trivial k-point set. e. unless you are looking at a molecule or very large system do not use the Gamma-point only.

3. Loop over a set of kinetic energy cut-off values.

   These should be a simple static calculation. Make sure that each of the calculations finishes successfully, otherwise you will not be able to compare results and check convergence.

4. Collect relevant data and check the convergence behavior.

Convergence of the kinetic energy cut-off for alpha Ce using the PBE functional and a 9x9x9 k-point grid.

In our example, we used a 9x9x9 k-point set. Looking at the example, we first of all see how smoothly the total energy varies with regard to the ENCUT parameter. In addition, it is important to note that VASP has a correction term (search for EATOM in the OUTCAR file) implemented which greatly improves the energy convergence (compare the black and red curves). Unfortunately, it also leads to non-variational convergence (i.e. the energy does not become strictly smaller with increasing cut-off) which may lead to some confusion. However, the correction term performs really well, and allows you to use a kinetic energy cut-off which is much lower than what you would need to use without. In this case, the default cut-off misses the reference energy by about 10 meV. Without the correction, a cut-off of about 540 eV (almost double) is needed. From ENCUT=300 to 800 eV, you observe a plateau, so using a higher cut-off will not improve the energy much. However, other properties, such as the calculated forces or the hessian may improve in this region. For these parameters a higher cut-off may be beneficial, and their convergence as function of ENCUT should be checked if important for your work.

2.     K-point set

Similar as for the kinetic energy cut-off, if you are working with a periodic system you should check the convergence of your k-point set. However, if you are working with molecules/clusters your Brillouin zone reduces to a single point, so your k-point set should only consist of the Gamma point and no convergence testing is needed. More importantly, if you use a larger k-point set for such systems (molecules/clusters) you introduce artificial interaction between the periodic copies which should be avoided at all cost.

For bulk materials a k-point convergence check has a similar setup as the basis set convergence check. The main difference being the fact that for these calculation the basis set is kept constant (VASP: ENCUT = default cut-off, manually set) and the k-point set is varied. As such, if you are new to quantum mechanical calculations, or start using a new code, you can combine the two convergence checks and study the convergence behavior on a 2D surface.

K-point convergence of alpha-Cerium using the PBE functional and ENCUT=500 eV.

In our example, ENCUT was set to 500 eV. It is clear that an extended k-point set is important for small systems, as the Gamma-point only energy can be off by several eV. This is even the case for some large systems like MOFs. An important thing to remember with regard to k-point convergence, is the fact that this convergence is not strictly declining, it may show significant oscillations overshooting and undershooting the converged value. A convergence of 1 meV or less for the entire system is a goal to aim for. An exception may be the most large systems, but even then one should keep in mind the size of the energy barriers in that system. Using flexible MOFs as an example which show a large-pore to narrow-pore transition barrier of 10-20 meV per formula unit, k-point convergence should be much below this. Otherwise your system may accidentally cross this barrier during relaxation.

The blue curve shows the number of k-points in the irreducible Brillouin zone. For standard density functional theory calculations (LDA and GGA, not hybrid functionals) this is a measure of the computational cost, as the k-points can be calculated fully independently in parallel (and yes the blue scale is a log-scale as well). Because the first orders of magnitude in accuracy are quickly crossed ( from Gamma to 6x6x6 the energy error goes from the order of eV to meV) while the number of k-points doesn’t grow that quickly (from 1 to 28). As a result, one often performs structure optimizations in a stepped fashion, starting with a coarse grid steadily increasing the grid (unless pathological behavior is expected… MOFs again…yes, they do leave you with nightmares in this regard).

3.     Conclusions

Convergence testing is necessary, in theory, for each and every new system you look into. Luckily, VASP behaves rather nicely, such that over time you will know what to expect and your convergence tests will reduce in size significantly and become more focused. In the examples above we used the total energy as a reference, but this is not always the most important aspect to consider. In some cases you should check the convergence as function of the accuracy of the forces. In that case you generally will end up with more stringent criteria as energy converges rather nicely and quickly.

May your convergence curves be smooth and quick.

Ball-and-stick representation of diamond in two different unit cells. Left: primitive unit cell containing two atoms. All atoms at the vertices are periodic copies of the same one. Right: Conventional cubic unit cell containing eight atoms. Atoms at opposing faces are periodic copies, while all atoms at the vertices are periodic copies of the same atom.

When performing electronic structure calculations on complex systems, you prefer to do this on systems with as few atoms as possible. Such periodic cells are called unit cells. There are, however, two types of unit cells: Primitive unit cells and conventional unit cells.

A primitive unit cell is the smallest possible periodic cell of a crystalline material, making it extremely suited for calculations. Unfortunately, it is not always the nicest unit cell to work with, as it may be difficult to recognize it’s symmetry (cf. the example of diamond on the right). The conventional unit cell on the other hand shows the symmetry more clearly, but is not (always) the smallest possible unit cell. To make matters complicated and confusing, people often refer to both types as simply “unit cell”, which is not wrong, but the term unit cell is for many uniquely associated with only one of the two types.

When you are performing calculations on diamond, the conventional cell isn’t that large that standard calculations become impossible, even on a personal laptop or desktop. On the other hand, when you are studying a Metal-Organic Framework like the UiO-66(Zr) which contains 456 atoms in its conventional unit cell, you will be very happy to use the primitive unit cell with ‘merely’ 114 atoms. Also the MIL-47/53 topology which generally is studied using a conventional unit cell containing 72/76 can be reduced to a smaller primitive unit cell of only 36/38 atoms. Just as for the diamond primitive unit cell, this MIL47/53 primitive unit cell is not a nice cubic cell. Instead you end up with a lattice having lattice angles of seventy-something degrees.

Reduction of the MIL-53 conventional cell to the primitive cell. The conventional cell is shown, extending slightly into the periodic copies. The primitive lattice vectors are shown as colored arrows. The folded primitive cell shows there was some symmetry-breaking in the hydroxy groups of the metal-oxide chain. Introducing some additional symmetry fixes this in the final primitive cell.

How to reduce a conventional unit cell to a primitive unit cell?

Before you start, and if you are using VASP, make sure you have the POSCAR file giving the atomic positions as Cartesian coordinates. (Using the HIVE-4 toolbox: Option TF, suboption 2 (Dir->Cart).)

If you do not use VASP, you can still make use of the scheme below.

1. Open your structure using VESTA, and save it as “VASP” file: POSCAR.vasp (File → Export Data → choose “VASP” as filetype, select Cartesian Coordinates (don’t select the Convert to Niggli reduced cell as this only works for perfect crystal symmetry)).
2. Open the file you just saved in a text editor (e.g. notepad or notepad++). The file format is quite straight forward.  The first line is a comment line while the second is a general scale-factor, which for our current purpose can be ignored. What is important to know is that the 3^rd, 4^th and 5^th lines give the lattice vectors (a, b, and c). The 6^th and 7^th line give the order, type and number of atoms for each atomic species (In VASP 5.x, the older VASP 4.x format does not have a 6^th line). The 8^th line should say “Cartesian”. From the 9^th line onward you get the atomic coordinates.
3. Choose 1 atom in your conventional cell which you are going to use as reference point.
4. Get the primitive unit cell lattice vectors by generating vectors from the reference atom. (cf. figure above) Using VESTA this can be done as follows:
   1. Open your conventional cell in VESTA (if you closed it after step 1).
   2. Use the distance selector (5th symbol from the top in the left-hand-side menu) and for each of the primitive lattice vectors, select the reference atom and it’s primitive copy.
      
   3. Subtract the “fractional coordinates” of the selected atoms provided by VESTA to get a “fractional” primitive vector (the primitive a vector will be called a^prim,frac )
   4. Multiply each of the conventional lattice vectors(a^conv, b^conv, and c^conv) with the corresponding component of the fractional primitive vector, and add the resulting vectors to obtain the new primitive vector:

      ⇒  a^prim = a_x^prim,frac a^conv + a_y^prim,frac b^conv + a_z^prim,frac c^conv

      So imagine that the lattice vectors of the MOF above are  a = ( 20, 0, 0),  b = ( 0, 15, 0), and c = ( 0, 0, 5). And the primitive fractional a vector is found to be a^prim,frac = ( 0.5, -0.5, 0.5). In this case the a^prim vector will become: a^prim = ( 10, 0, 0 ) + ( 0, -7.5, 0 ) + ( 0, 0, 2.5) = (10, -7.5, 2.5).

5. Replace the conventional lattice vectors in the POSCAR.vasp file (cf. step 2) with the new primitive lattice vectors. Save the file.
6. Open the POSCAR.vasp in VESTA. If everything went well, and the conventional cell wasn’t the real primitive cell already, you should see a nice new primitive cell with the equivalent atoms perfectly overlapping one-another. This is also the reason to have your starting geometry in Cartesian coordinates. If you would have your atomic positions as fractional coordinates this first check will not work at all. Furthermore, you would need to calculate the new fractional coordinates of the atoms in the primitive unit cell. If all is well, you can close POSCAR.VASP in VESTA. (If something is wrong: either you did something wrong, and you should start again, or it wasn’t actually a super cell of a primitive cell you started to construct.)
7. Get the atoms of the primitive unit cell.
   1. Because our atomic positions are in Cartesian coordinates in our initial geometry file, we now just need to make a list of single copies of equivalent atoms. Using VESTA (the original structure file you still have open from step 1) you can click on each atom you wish to keep and write down their index (this is the first number you find on the line with Cartesian coordinates) … For example: In case of the MIL53-MOF you can select all metal and oxygen atoms of 1 chain, and two linker molecules.
   2. Remove all superfluous atoms (i.e. those of which you didn’t write down the index) from the POSCAR.vasp (using your text-editor). You may want to make a backup of this file before you start :-).
   3. Update the number of atoms on the 7^th line of the POSCAR.vasp file, and check that the number is correct. The conventional cell should have had an integer multiple of the number of atoms in the primitive cell. Save the final structure as POSCAR_final.vasp .
   4. The POSCAR_final.vasp should contain both the new lattice vectors and a list of atoms for a single primitive unit cell. Check this by opening the file using VESTA, and make sure you didn’t remove too many or too few atoms. If not, go back to step (a) and double check. (If you are using the HIVE4-toolbox you can first transform POSCAR_final.vasp back to direct coordinates as this may make atoms visible which nicely overlap in Cartesian coordinates: Option TF, suboption 1 (Cart->Dir)  )
8. Congratulations you have constructed a primitive cell from a conventional cell.

 As you can see, the method is quite simple and straight forward, albeit a bit tedious if you need to do this many times.

Enjoy your primitive unit cell!

PS: Small remark for those new to VESTA. You can use delete atoms in VESTA and store your structure again. This is useful if you want to play with a molecule. Unfortunately for a solid you need also to get new lattice vectors, which did not happen. As a result you end up with some atoms floating around in a periodically repeated box with the original lattice parameters. Steps 1-5 given above provide a simple way of not ending up in this situation, but require some typing on your part.

PS 2: The opposite transformation, from a primitive unit cell to a conventional unit cell, using VESTA, is shown in this youtube video.