Quantum Espresso to Perturbo

Before running electron dynamics calculations using pertubo.x, the user needs to carry out electronic and phonon calculations, with DFT and DFPT respectively. At present, Perturbo can read the output of DFT and DFPT calculations done with Quantum Espresso (QE). Once the relevant output files have been obtained from QE and Wannier90 (W90), the first step is to use qe2pert.x to compute the e-ph matrix elements on a coarse \(\mathbf{k}\) and \(\mathbf{q}\) point Brillouin zone grid, to obtain e-ph matrix elements in Wannier function basis, and to store the data into the HDF5 format for perturbo.x to read. The generation of this HDF5 file, called ‘prefix’_epr.h5, is discussed in this section of the manual.

The preparation stage consists of five steps:

1. Run a self-consistent (SCF) DFT calculation
2. Run a phonon calculation using DFPT
3. Run a non-SCF (nSCF) DFT calculation
4. Run Wannier90 to obtain Wannier functions
5. Run qe2pert.x

In the following, we use silicon as an example. The input files for QE and W90 are in the directory “example02-silicon-qe2pert/pw-ph-wann” of the repository of perturbo-examples-light. As a reference, we also provide the results in a directory called “References”.

Step 1: SCF calculation

Here is the scf.in QE input file:

&CONTROL
  prefix = 'si'
  calculation = 'scf'
  wf_collect = .true.
  outdir='./tmp'
  pseudo_dir='../pseudo'
/
&SYSTEM
  ibrav = 2
  celldm(1) = 10.264
  nat = 2
  ntyp = 1
  ecutwfc = 40.0
/
&ELECTRONS
  conv_thr = 1.0d-15
  mixing_mode = 'plain'
  mixing_beta = 0.7
  diagonalization = 'david'
  diago_full_acc = .true.
/
&IONS
/
&CELL
press_conv_thr = 0.01
/
ATOMIC_SPECIES
  Si   28.085  Si_DOJO_PBEsol.upf
ATOMIC_POSITIONS crystal
Si  0.00000000  0.00000000  0.00000000
Si -0.25000000  0.75000000 -0.25000000
K_POINTS (automatic)
 20 20 20 0 0 0

Run an SCF calculation (pw.x QE executable) and obtain the QE ‘prefix’.save directory. In this case, we obtain ./tmp/si.save, which is needed for phonon and nSCF calculations.

Step 2: phonon calculation

Here is the ph.in input file:

Phonons on a uniform grid
&inputph
  verbosity='debug'
  tr2_ph=1.0d-17
  prefix='si'
  ldisp=.true.
  epsil=.true.
  lqdir=.true.
  outdir='./tmp'
  fildyn  = 'si.dyn.xml'
  fildvscf = 'dvscf'
  nq1=8, nq2=8, nq3=8,
/

Copy the QE ‘prefix’.save directory from the SCF run to the current directory:

$ cp -r ../scf/tmp ./

In this step, make sure that the number of \(\mathbf{q}\) points is commensurate with the number of \(\mathbf{k}\) points used in the nSCF and Wannierization calculations. For example, a \(\mathbf{q}\) grid of 8x8x8 can be used with a wannierization \(\mathbf{k}\) grid of 8x8x8 or 16x16x16, but not with a 10x10x10 or 12x12x12 grid.

Run the PHonon QE calculation (ph.x QE executable).

Once the calculations are done, we collect the phonon data into a directory called save, created by running the shell script ph-collect.sh.

$ ./ph-collect.sh

or

$ bash ph-collect.sh

The save directory contains all the information needed for PERTURBO to interface with QE. It includes the dynamical matrix files, phonon perturbation potentials, and the patterns.

Since phonon calculations using DFPT can be computationally expensive, it is often useful to estimate the number of irreducible \(\mathbf{q}\) points before running the phonon calculation (run a Gamma-point phonon calculation). Note however that this step is optional.

Step 3: nSCF calculation

We now run the nSCF calculations needed to generate the wavefunctions on the full \(\mathbf{k}\) point grid, which we’ll need both for generating Wannier functions with Wannier90 and for forming the coarse-grid e-ph matrix elements in Perturbo.

Here is the nscf.in input file (truncated):

&CONTROL
  prefix = 'si'
  calculation = 'nscf'
  wf_collect = .true.
  outdir='./tmp'
  pseudo_dir='../pseudo'
/
&SYSTEM
  ibrav = 2
  celldm(1) = 10.264
  nat = 2
  nbnd=16
  ntyp = 1
  ecutwfc = 40.0
/
&ELECTRONS
  conv_thr = 3.0d-12
  mixing_mode = 'plain'
  mixing_beta = 0.7
  diagonalization = 'david'
  diago_full_acc = .true.
/
&IONS
/
&CELL
press_conv_thr = 0.01
/
ATOMIC_SPECIES
  Si   28.085  Si_DOJO_PBEsol.upf
ATOMIC_POSITIONS crystal
Si  0.00000000  0.00000000  0.00000000
Si -0.25000000  0.75000000 -0.25000000
K_POINTS crystal
512
  0.00000000  0.00000000  0.00000000  1.953125e-03
  ...

Make sure that the number of k points is commensurate with the number of \(\mathbf{q}\) points used for phonons, otherwise, qe2pert.x will stop. Remember to copy the QE ‘prefix’.save directory from the SCF calculation the current directory:

$ cp -r ../scf/tmp ./

The k points lists could be produced by using kmesh.pl in ‘utility’ folder of Wannier90:

$ kmesh.pl 8 8 8 >> nscf.in

Then run the nSCF calculation with QE (pw.x executable).

Step 4: Wannier90 calculation

The directory contains two input files, one for wannier.x and the other for pw2wannier90.x.

Here is the pw2wan.in input file:

&inputpp
  outdir='./tmp'
  prefix = 'si'
  seedname = 'si'
  spin_component = 'none'
  write_mmn = .true.
  write_amn = .true.
  write_unk = .false.
/

And here is the si.win input file (truncated):

begin projections
 Si:sp3
end projections
guiding_centres=true

num_bands = 16
num_wann = 8

iprint = 2
dis_num_iter =  500
dis_win_min =  -100.000
dis_win_max =   17.200
dis_froz_min = -100.000
dis_froz_max  = 9.000
num_iter  =   10000
mp_grid : 8 8 8

begin unit_cell_cart
bohr
-5.1320  0.0000  5.1320
 0.0000  5.1320  5.1320
-5.1320  5.1320  0.0000
end unit_cell_cart

write_u_matrices = .true.
write_xyz = .true.

#restart = plot
BANDS_PLOT = TRUE
BANDS_PLOT_FORMAT = gnuplot
BANDS_NUM_POINTS = 100

BEGIN KPOINT_PATH
L 0.500 0.500 0.500  G 0.000 0.000 0.000
G 0.000 0.000 0.000  X 0.500 0.000 0.500
X 0.500 0.000 0.500  W 0.500 0.250 0.750
W 0.500 0.250 0.750  K 0.375 0.375 0.750
K 0.375 0.375 0.750  G 0.000 0.000 0.000
END KPOINT_PATH

begin atoms_frac
 Si    0.00000   0.00000   0.00000
 Si   -0.25000   0.75000  -0.25000
end atoms_frac

begin kpoints
 0.0000000  0.0000000  0.0000000
 ...
end kpoints

In the input file si.win, we instruct Wannier90 to write two important quantities for qe2pert.x, the \(U(\mathbf{k})\), \(U^{\text{dis}}(\mathbf{k})\) matrices and the position of the Wannier function centers, using: write_u_matrices=true and write_xyz=true.

We create tmp directory:

$ mkdir tmp

and change into it. We soft link to the QE ‘prefix’.save directory obtained in the nSCF calculation:

$ cd tmp
$ ln -sf ../../nscf/tmp/si.save 

We then run Wannier90. This step consists of three runs:

• 1) Run wannier90.x -pp si to generate a list of the overlaps (the si.nnkp file, more details here).

• 2) Run pw2wannier90.x < pw2wan.in.

• 3) Run wannier90.x si (the si.win input file will be implicitly used).

The important output files for qe2pert.x are si_u.mat, si_u_dis.mat, and si_centres.xyz. For disentangled bands, there would be no ‘prefix’_u_dis.mat. We encourage the user to become familiar with Wannier90 to run this step for different materials.

The user has to run Wannier90 3.0 or higher, since otherwise the \(U\) matrices will be printed out.

Step 5: Running qe2pert.x

We are ready to compute the e-ph matrix elements on the coarse \(\mathbf{k}\) point (determined by the nSCF step) and \(\mathbf{q}\) point (determined by the phonon step) Brillouin zone grids. First, copy or link the electronic and phonon calculation results to the current directory.

$ cd qe2pert
$ mkdir tmp
$ cd tmp 

$ #link to the nscf .save directory
$ ln -sf ../../pw-ph-wann/nscf/tmp/si.save
$ cd ../

$ #link to the wannier information
$ ln -sf ../pw-ph-wann/wann/si_u.mat
$ ln -sf ../pw-ph-wann/wann/si_u_dis.mat
$ ln -sf ../pw-ph-wann/wann/si_centres.xyz

Here we show the input file (qe2pert.in) for the executable qe2pert.x:

&qe2pert
 prefix='si'
 outdir='./tmp'
 phdir='../pw-ph-wann/phonon/save'
 nk1=8, nk2=8, nk3=8
 dft_band_min = 1
 dft_band_max = 16
 num_wann = 8
 lwannier=.true.
 load_ephmat = .false.
 system_2d = .false.
/ 

The description of the input parameters:

• prefix: needs to be the same as the prefix used in the input files for QE.
• outdir: contains the save directroy obtained from the nSCF calculations. The calculated e-ph matrix elements will be stored in this directory.
• phdir: is the save directory inside which we collected all the phonon information.
• nk1, nk2, nk3: are the number of \(\mathbf{k}\) points along each direction used in the nSCF and Wannier90 calculations.
• dft_band_min and dft_band_max: determine the range of bands we are interested in, and should be the same as the values used in the Wannierization process. For example, if we used 40 bands in the nSCF calculation and we excluded bands 1-4 and 31-40 in the Wannierization, then dft_band_min=5 and dft_band_max=30.
• num_wann: the number of Wannier functions.
• lwannier: a logical flag. When it is .true., the e-ph matrix elements are computed using the Bloch wave functions rotated with the Wannier unitary matrix; if .false., the e-ph matrix elements are computed using the Bloch wave functions, and the e-ph matrix elements are then rotated using the Wannier unitary matrix. By default, it is .true. to reduce computational cost.
• load_ephmat: a logical flag. If .true., reuse e-ph matrix elements in Bloch function basis computed previously. This is useful if you want to test different Wannier bases. For example, you could first run qe2pert.x with lwannier=.false., and then rerun qe2pert.x with lwannier=.false. and load_ephmat=.true. with different Wannier unitary matrix.
• system_2d: if the materials is two-dimensional, so that in one direction only one \(\mathbf{k}\) point is used, set it to .true.; the default is .false..

Now we are ready to run the e-ph matrix elements:

export OMP_NUM_THREADS=4
$ mpirun -np 2 qe2pert.x -npools 2 -i qe2pert.in > qe2pert.out

This task is usually time-consuming on a single core, but it can be made much faster (minutes) on multiple cores. The executables qe2pert.x employ hybrid parallelization (MPI plus OpenMP), e.g. 2 MPI processes and each process span 4 OpenMP threads in this example.

To speed up the calculations, the users could increase the number of OpenMP threads and MPI processes. Threads with OpenMP are particularly useful when the RAM (memory) of computing nodes is limited. The memory consumption reduces to minimum when using 1 MPI process per node and setting OMP_NUM_THREADS to the number of cores per node.

Once the calcalculation has completed, we obtain the output file si_epr.h5, which is an HDF5 database with all the information needed to run perturbo.x (which is described in the next section).

Structure of epr.h5 HDF5 file

The HDF5 file ‘prefix’_epr.h5 has the following data structure:

The ‘prefix’_epr.h5 file structure can be schematically represented as follows: