Notes on Quantum ESPRESSO

Here I collect some notes on how to use Quantum ESPRESSO version 5.** to perform some standard tasks.
There is no guarantee that this will work as you desire, so use them at your own risk.

Link to the official site: Quantum ESPRESSO

Index

Calculate bands
Calculate DOS
Calculate epsilon
Calculate PDOS
Plot wavefunction

Calculate Bands

First one has to run a SCF calculation.
Then the type of calculation has to be changed to 'bands', and a set of k-points has to be given at the end:

....
calculation = 'bands' ,
outdir = './out/' ,
....
occupations = 'smearing' ,
....
K_POINTS crystal
119
0.333333333 0.333333333 0.000000000 1.000000000
0.323333333 0.323333333 0.000000000 1.000000000
0.313333333 0.313333333 0.000000000 1.000000000
0.303333333 0.303333333 0.000000000 1.000000000
....

A program has to be run in order to separate the band values from the rest of the data.
So an input file has to be created, bands.in

Example of file "bands.in":
&BANDS
prefix = 'mt_gipaw',
outdir='./out/',
filband = 'bands.dat',
/
and the program must be run (several processors):
mpirun -np 4 bands.x < bands.in > bands.out
Single processor:
bands.x < bands.in > bands.out

It has to run on the same number of processors as the first run.
To avoid this, read instructions of parallel running and add option to put all files in the same place.

A file bands.dat is generated, with the k-points and the energies, but still not very readable.
So another program has to be run to transform this to a readable thing.

plotband.x bands.dat

It asks several things:

Range: -19.1060 7.1480eV Emin, Emax > -20 8
output file (xmgr) > bandas.xmgr
output file (ps) > bandas.ps
Efermi > -4
deltaE, reference E (for tics) 5 -4

It generates a ps picture and the range and deltaE and reference is for that picture.

Calculate DOS

After calculating the scf for the material (with occupations='tetrahedra'), one has to run a nscf calculation where we have to specify:

For a DOS calculation, you should specifya uniform grid of points. For DOS calculations you should choose occupations='tetrahedra', together with an automatically generated uniform k-point grid (card K POINTS with option 'automatic'). ....
occupations = 'tetrahedra' ,
....
and the number of k-points for the sampling of DOS:
....
K_POINTS automatic
48 48 48 0 0 0

EXAMPLE:
&CONTROL ....
calculation = 'nscf' ,
restart_mode = 'from_scratch' ,
outdir = './out/' ,
....
occupations = 'tetrahedra' ,
....
K_POINTS automatic
48 48 1 0 0 0

Then one has to run the program dos.x < dos.in (single processor)
EXAMPLE of dos.in:

&dos
prefix = 'gipaw' ,
outdir = './out/' ,
fildos='dos.dat',
Emin=-10.0, Emax=5.0, DeltaE=0.01
/

DeltaE is the energy stepping, between Emin and Emax.
WARNING: degauss is gaussian broadening, Ry (not eV!) read from the input data file (they will be the same used in the pw.x calculations)
WARNING: occupations = 'tetrahedra' ,

The total DOS (states/eV plotted vs E in eV) is written:
states/eV/(unit cell)

Usually:
DOS(E) dE = number of energy levels in the energy range from E and E+dE and according to this definition
\int_E0^E1 DOS(E) dE = total number of states between E0 and E1 (adimensional).
This is what the dos.x executable included in Quantum-ESPRESSO computes. According to the above definition:
DOS(E) = \sum_n \int delta(E - E_n(k_x,k_y,k_z)) dk_x dk_y dk_z *V / (4 \pi^3)

Calculate epsilon

Instructions for producing data on dielectric constant

Make a normal scf calculation.
Then run a nscf calculation where you include:
nosym = .true. ,
noinv = .true. ,

(This way we can use automatic k-point generation)
K_POINTS automatic
50 50 50 0 0 0

Then you run the program epsilon.x (or epsFR.x, for noncolinear calculations), with the same number of processors:

nohup mpirun -np 2 epsilon.x < eps.in > eps.out &

Example of file eps.in

&inputpp
outdir='./out/',
prefix='MoS2',
calculation='eps',
/
&energy_grid
smeartype='gauss',
intersmear=0.136d0,
intrasmear=0.0d0,
wmax=10.0d0,
wmin=0.0d0,
nw=500,
shift=0.0d0,
/

WARNING: In this file, energy units are eV.

intersmear is the broadening for interband transitions. If it is too small, a number of spurious peaks due to numerical reasons will show up. I don't know any rule how to choose it. What should be done is try several broadenings and see which features remain and which move. A value around 0.1 eV sounds reasonable.

intrasmear is the broadening for intraband transitions. It is nedded when you have a metal, and it should have a value if you use smear in the scf calculations. Again, a value around 0.1 eV sounds reasonable.

In particular systems, like graphene, these calculations have to be done very carefully.
Since the fermi surface is reduced to a dot, an extremelly high density of k-points have to be used to achieve accuracy near the fermi level (that is, for low transition energies). 750x750 is a good number, but with 520x520 it is possible to get a good result.
It is mandatory to use intrasmear; a good value is between 0.15 and 0.2 eV.

If the system is bi-dimensional, one has to multiply the value by the height of the supercell (in SI), to get the right value.

Calculate PDOS

Perform scf calculation.
Perform nscf calculation:
if using full relativistic with spin orbit, you have to do a nscf calculation with the options:
nosym = .true. ,
occupations = 'smearing' ,
K_POINTS automatic
16 16 1 0 0 0
with shift = 0 0 0 !
* The tetrahedron method is presently not implemented.

Perform a pdos calculation:
Can be with single processor.
run: projwfc.x < pdos.in > pdos.out
example of pdos.in:

&PROJWFC
prefix = 'MX2' ,
outdir = './out/' ,
filpdos = 'Zr-pdos' ,
degauss = 0.01 ,
DeltaE = 0.01 ,
Emin = -5 ,
Emax = 5 ,
/

Make sure the Ef is between Emin and Emax!
In the non-collinear, spin-orbit case (i.e. if there is at least one fully relativistic pseudopotential) wavefunctions are projected onto eigenstates of the total angular-momentum.
Projected DOS are written to file
"filpdos".pdos_atm#N(X)_wfc#M(l_j),
where N = atom number , X = atom symbol, M = wfc number, l=s,p,d,f and j is the value of the total angular momentum.

Example: Zr-pdos.pdos_atm#3(Zr)_wfc#2(p_j1.5)

In this case the format is:
E LDOS(E) PDOS_1(E) ... PDOS_2j+1(E)
...
where LDOS = \sum m=1,2j+1 PDOS_m(E)
and PDOS_m(E) = projected DOS on atomic wfc with component m

Example:
# E(eV) ldos(E) pdos(E)_1 pdos(E)_2 pdos(E)_3 pdos(E)_4
-5.000 0.480E-01 0.872E-02 0.153E-01 0.153E-01 0.872E-02
-4.990 0.483E-01 0.862E-02 0.155E-01 0.155E-01 0.862E-02
-4.980 0.483E-01 0.848E-02 0.157E-01 0.157E-01 0.848E-02

All DOS(E) are in states/eV plotted vs E in eV

::: Orbital Order
Order of m-components for each l in the output:
1, cos(phi), sin(phi), cos(2*phi), sin(2*phi), .., cos(l*phi), sin(l*phi)
where phi is the polar angle:x=r cos(theta)cos(phi), y=r cos(theta)sin(phi)
This is determined in file flib/ylmr2.f90 that calculates spherical harmonics.

for l=1:
1 pz (m=0)
2 px (real combination of m=+/-1 with cosine)
3 py (real combination of m=+/-1 with sine)

for l=2:
1 dz2 (m=0)
2 dzx (real combination of m=+/-1 with cosine)
3 dzy (real combination of m=+/-1 with sine)
4 dx2-y2 (real combination of m=+/-2 with cosine)
5 dxy (real combination of m=+/-1 with sine)

Plot wavefunction

Instructions how to plot the wavefunctions using xcrysden

First run a scf calculation.
Then run a nscf calculation with just the k-points you are interested.
Then, for each k-point and each band, you have to run the pp.x code.

The pp input file can be divided in two operations:
1) extracting the information to a file
2) extracting from that file the more specific information wanted, and save it to a formated file

Example of pp input file for the first operation:
&INPUTPP
prefix = 'MoS2' ,
outdir = './out' ,
filplot = 'M18-0' ,
plot_num = 7,
spin_component = 0 ,
kpoint = 4,
kband = 18,
/

filplot is the file where the information will be saved.
plot_num is the type of plot, and is related to the spin_component
kpoint is the k-point where the calculation is made; the number indicates the order number of the nscf calculation
kband is the number of the band, starting from the lowest energy band.

Example of pp input filefor the second operation:

&INPUTPP
/
&PLOT
nfile = 1 ,
filepp(1) = './M18-0',
weight(1) = 1.0,
fileout = 'M18-0-plot' ,
iflag = 3 ,
output_format = 5 ,
/

These two operations can be put in the same file.

The output file M18-0-plot can be now read by xcrysden, using file-Open Structure-open XSF structure.
Then proceed via the Tools-Data Grid menu to show the image.
Find the ok button to proceed (you may have to increase the size of the window to see it).
Then a window with parameters that can be controlled, including switching on the image.