Fatbands with spin-orbit coupling

A recent question touched on the case of fatband representation when spin-orbit coupling is included in the calculation (lspinorb = .true.). I assumed that this would be straightforward and indeed there is no problem up to the last step, where plotband.x should read the band energies and projections. It seems that in QE 6.5 the projection file is not properly read

 Input file > bands.dat
Reading 30 bands at 401 k-points
file with projections not compatible with bands
Range: -5.6460 34.2420eV Emin, Emax >

However the projection file is correctly written as we can identify the (j,mj) blocks such as

1    1  Pb  5D     1    2  2.5  -2.5

To circumvent this incompatibility of plotband.x one can simply read the energy eigenvalues from the bands.x outfile and relative weight of the wave-function from the projection file, where every block looks something likes this:

    9    1  Pb  5D     2    2  1.5   0.5
       1       1        0.5475316611
       1       2        0.2162196288
       1       3        0.1468317944
       1       4        0.0885236565
       1       5        0.0000000000
       1       6        0.0000000000
     401      29        0.0000000000
     401      30        0.0000000000

The first column is the k-point index, followed by the band and the weight for each mj state. If we want to plot the fatbands for all the d-states we have to identify the “block” number (i.e. the preceding 9 in the example above), collect all the d-states projections and plot them on top of the d-bands. I found a script in the usergroup and the origin seems to trace back to https://yyyu200.github.io/DFTbook/ but unfortunately I can’t read Chinese and cannot verify that information.

And there we have it: fatband representation for the d-bands of bulk lead with spin-orbit coupling.

A script that does all of this automatically is as always included, the pseudopotential used is the fully relativistic version of the ONCV pseudopotential which seems to be giving reliable results according to SSSP, however I did not converge the calculation as this seemed precise enough for a quick demonstration.

I hope that helps!

Five easy ways to visualize xsf volume grids

Post-processing of the total charge density allows the visualization of many quantities such as charge density, spin-polarization, electron localization and many more (see the option plot_num in the manual). In many cases the resulting data will be a 3D grid representation (voxel grid) in either the CUBE format or XSF – I prefer the latter because I am more familiar with using it.  

For this example I generated an xsf file of the spin-polarization (plot_num=6) of a molecule (I used nickelocene as an example as I had it readily at hand). By default, pp.x will generate the xsf on the entire FFT grid, which might be bigger than necessary. To avoid this, you can use “bspline” interpolation together with output_format = 3 which seems to work without any issue in QE 6.5. 

The format of the xsf file is pretty self-explanatory and also detailed here. The header contains information about the cell dimensions, the coordinates of the atoms in the cell, followed by the grid specifications and the grid data “one pixel at a time”, which are written as explained in the link as follows:

write(*,*) $ (((value(ix,iy,iz),ix=1,nx),iy=1,ny),iz=1,nz)

Now, let’s move on to how to visualize this data. The most “natural” way would probably be using xcrysden, which is excellent when dealing with xsf files, but I wanted to use MATLAB first to get a better feeling of how to read the xsf raw data so let’s look at a simple MATLAB script to read and visualize xsf voxel data first. 

1) MATLAB (https://www.mathworks.com/

Reading in the voxel grid is rather easy; I opted to look for the keyword DATAGRID_3D_UNKNOWN first and then simply read in all following lines starting 6 lines after the beginning of the datagrid. There might be much better ways of doing this, but it was all I could come up with in 5 minutes or less. In the end we have to reshape the data which can be easily achieved by MATLAB’s own reshape routine. Isosurfaces can be constructed (and simultaneously plotted) by the isosurface() function. The script to read the xsf file and create the following plot is attached at the end of this post.

I have centered the molecule at (0/0/0) for plotting purposes and I would say the result isn’t bad for such a simple approach! If someone wants to manipulate or use the data to calculate volume overlaps or similar then MATLAB is definitely convenient, but when it comes to creating pretty plots I am not sure how much this can be improved. So let’s move to option  

2) XCRYSDEN (http://www.xcrysden.org/

It does not get much easier to read and visualize xsf files than with xcrysden. Simply call xcrysden with the — xsf  switch (I am using xcrysden 1.6.2 on Ubuntu), press a to activate anti-aliasing to make it  a bit prettier and click on Tools -> Data Grid -> OK; input the isovalue, check to render positive and negative isovalues and you are ready to go; Unfortunately exporting the resulting file to an image is a bit limited and the pixel based images (png, TIFF etc.) are not the prettiest and their resolution is limited to the screen resolution. However, this is definitely a fast and convenient way. 

3) VESTA (http://jp-minerals.org/vesta/en/

Slightly nicer plots can be made with next to no effort in VESTA as well. Just drag and drop the xsf file into a recent version of VESTA (I use version 3.4.4 dated March 2018) and adjust the settings under “Properties -> Isosurfaces” to your liking. What I particularly like in VESTA is the easy slicing using the “Boundary -> Cutoff planes -> new) function. To get a better resolution, export the raster image (File -> Export raster image) at larger sizes, then shrink it down.

4) OVITO (https://www.ovito.org/

I have only used this one for a short period of time compared to the others but OVITO (Open Visualization Tool) is hands down my favorite. The following plots were made with version 3.0.0 and there is also a “Pro” version available that will come with license costs in the near future. Simply drag-and-drop the xsf file into the main window, then add a “Create isosurface” modifier (Add modification -> Visualization -> Create Isosurface); do this twice, once for positive and once for negative values of the isosurface and you are ready to export the raster image with a render engine of your choice.  Here I have used Tachyon which already gives quite nice results.  

The GUI offers a huge variety of customizations for atom color, radius, bond customization, slicing (Add Modification -> Modification -> Slice) making this tool really fun to work with. But if you want to go one step further you can use… 


OK, I admit, this is not really another tool but rather the python scripting interface to OVITO. Nevertheless, this is an amazing addition that makes repeated tasks much less annoying. When using the script interface you can further take advance of numpy which makes “doing math” on the xsf file very easy. This approach really shines when you have to create a lot of similar plots for a publication and you want to make sure that plotting is consistent across these plots (for example atomic radii or colors). I have attached the python script used to create this plot which I hope is basically self-explanatory. 

 If you have any questions drop them in the comments and I will try to help you out! 

Happy Plotting! 

Work function calculation of MgO in QE 6.2.1

Using a slab-vacuum system and the dipole correction as previously described we aim to calculate the work-function (WF) of MgO. The experimental value for the WF of a low-index surface such as (200) is about 5 eV according to Choi et al.


By constructing a 6 layer MgO slab padded by 6 layers of vacuum we create an artificial surface that is infinitely repeated by periodic boundary conditions. The dipole correction is applied by setting tefield and dipfield to true. We know that the position of the dipole is critical as the discontinuity has to fall into the vacuum region. By moving the dipole across the vacuum region it can be shown that if the dipole is too close to the slab the dipole correction will fail.


If placed correctly the electrostatic potential will be flat in the vacuum region. Two failing cases are also shown by dashed lines. The slab is repeated to the negative region to emphasize the region in which the dipole correction is physical which extends from emaxpos+eopreg-1 (in units of entire cell extending from [0-1]). The difference between vacuum level (here set to 0) and the flat region of the electrostatic potential is about 5 eV which is in acceptable agreement with experiment. Find the script to run this job here.