1

ZnS (zincblende structure) electronic structure calculation

using Quantum Espresso

Before you begin a project like this, identify the research question and compound(s) of interest.
Review what is known (experimentally) about the compound or the class of compounds to which it
belongs. If it has been studied before, what level of theory was brought to bear in the analysis.
What unique tools are available to you to build on this foundation? Are you/your lab the rigth group
to take on this project or is there a willing collaborator in the wings able - with resources and sta -
to help?

In this example, I looked at a recent paper examining the electronic structure of doped zinc
sul de.1

1. After identifying the material of intest, obtain - or generate - crystallographic coordinates for the
structure. The simplest way to do this is obtain a .cif from the CCSD or ICSD. Again, in this
case, I obtained the .cif for ZnS (zincblende).2

2. Generate PWscf.in le for pw.x calculation using either cif2cell (and a suitable .cif) and then
modifying the resulting output le or the Quantum ESPRESSO input generator and structure
visualizer (https://www.materialscloud.org/work/tools/qeinputgenerator) available at
materialscloud.org. In addition to a correctly formatted input le, the QE input generator takes
a .cif and arrives at informed parameter settings not obvious to the novice. Moreover, the
PWscf.in includes pseudopotential les in a .zip folder for easy transfer to your workstation,
pegasus, actinide or external resource. If you save this le, you can scp it to another computer
on the network.

Note that materialscloud.org developed novel methods for the calculation that you are about to run
and should be cited accordingly.3

1 D’Amico, P., Calzolari, A., Ruini, A. et al. New energy with ZnS: novel applications for a standard
transparent compound. Sci. Rep. 2017 7, 16805. https://doi.org/10.1038/s41598-017-17156-w
2E.A. Jumpertz, Zeitchrift fur Elektrochemie 1955, 59, 419. https://www.crystallography.net/cod//
1100043.html, ICSD 60378
3 If you use the results of this tool in a publication, please cite the following works:
• SSSP (for the pseudopotential library) G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet and
N. Marzari, npj Computational Materials 4, 72 (2018). WEB: http://materialscloud.org/sssp.
• Pseudopotentials:
◦ s_pbesol_v1.4.uspp.F.UPF,
◦ zn_pbesol_v1.uspp.F.UPF,
◦ from GBRV: K. F. Garrity, J. W. Bennett, K. M. Rabe, and D. Vanderbilt, Comput.
Mater. Sci. 81, 446 (2014).
◦ DOI: 10.1016/j.commatsci.2013.08.053, WEB: http://www.physics.rutgers.edu/gbrv,
LICENSE: GNU General Public License (version 3).

Revised 7/10/24

fi
fi
fi
fi
fi
fi
ff
 1

&CONTROL
calculation = 'scf'
etot_conv_thr = 8.0000000000d-05
forc_conv_thr = 1.0000000000d-04
outdir = './out/'
prefix = 'aiida'
pseudo_dir = './pseudo/'
tprnfor = .true.
tstress = .true.
verbosity = 'high'
/
&SYSTEM
degauss = 1.4699723600d-02
ecutrho = 7.2000000000d+02
ecutwfc = 9.0000000000d+01
ibrav = 0
nat = 8
nosym = .false.
ntyp = 2
occupations = 'smearing'
smearing = 'cold'
/
&ELECTRONS
conv_thr = 1.6000000000d-09
electron_maxstep = 80
mixing_beta = 4.0000000000d-01
/
ATOMIC_SPECIES
S 32.065 s_pbesol_v1.4.uspp.F.UPF
Zn 65.38 zn_pbesol_v1.uspp.F.UPF
ATOMIC_POSITIONS crystal
Zn 0.0000000000 0.0000000000 0.0000000000
Zn 0.0000000000 0.5000000000 0.5000000000
Zn 0.5000000000 0.0000000000 0.5000000000
Zn 0.5000000000 0.5000000000 0.0000000000
S 0.2500000000 0.2500000000 0.2500000000
S 0.2500000000 0.7500000000 0.7500000000
S 0.7500000000 0.2500000000 0.7500000000
S 0.7500000000 0.7500000000 0.2500000000
K_POINTS automatic
6 6 6 0 0 0
CELL_PARAMETERS angstrom
5.4145000000 0.0000000000 0.0000000000
0.0000000000 5.4145000000 0.0000000000
0.0000000000 0.0000000000 5.4145000000

3. Run pw.x to perform an SCF calculation and generate the Kohn-Sham states (DFT solutions to
Schrödinger Wave Equation). For big molecules, this should be done in parallel using mpirun:

/opt/intel/oneapi/mpi/2021.1.1/bin/mpirun -np 16 /opt/qe-7.1/bin/pw.x -i

pwscf.in > pwscf.out &

 2

At the end of pwscf.out you will nd the following:

the Fermi energy is 7.6352 ev

! total energy = -1941.49611650 Ry

estimated scf accuracy < 8.0E-12 Ry
smearing contrib. (-TS) = 0.00000000 Ry
internal energy E=F+TS = -1941.49611650 Ry

The total energy is F=E-TS. E is the sum of the following terms:

one-electron contribution = -1821.23804867 Ry
hartree contribution = 936.94398012 Ry
xc contribution = -200.48540163 Ry
ewald contribution = -856.71664632 Ry

convergence has been achieved in 12 iterations

. . .

This run was terminated on: 10:38:27 9Jul2024

=------------------------------------------------------------------------------=
JOB DONE.
=------------------------------------------------------------------------------=

4. The pwscf.in is used to generate pw.nscf.in and that le modi ed (change 'scf' to 'nscf',
increase the number of k-points) for a non-SCF calculation wherein a more dense k-mesh is
used. This produces a lower and more accurate energy.

&CONTROL
calculation = 'nscf'
etot_conv_thr = 8.0000000000d-05
forc_conv_thr = 1.0000000000d-04
outdir = './out/'
prefix = 'aiida'
pseudo_dir = './pseudo/'
tprnfor = .true.
tstress = .true.
verbosity = 'high'
/
&SYSTEM
degauss = 1.4699723600d-02
ecutrho = 7.2000000000d+02
ecutwfc = 9.0000000000d+01
ibrav = 0
nat = 8
nosym = .false.
ntyp = 2
occupations = 'smearing'
smearing = 'cold'
/
&ELECTRONS
conv_thr = 1.6000000000d-09
electron_maxstep = 80

fi
fi
fi
 3

mixing_beta = 4.0000000000d-01
/
ATOMIC_SPECIES
S 32.065 s_pbesol_v1.4.uspp.F.UPF
Zn 65.38 zn_pbesol_v1.uspp.F.UPF
ATOMIC_POSITIONS crystal
Zn 0.0000000000 0.0000000000 0.0000000000
Zn 0.0000000000 0.5000000000 0.5000000000
Zn 0.5000000000 0.0000000000 0.5000000000
Zn 0.5000000000 0.5000000000 0.0000000000
S 0.2500000000 0.2500000000 0.2500000000
S 0.2500000000 0.7500000000 0.7500000000
S 0.7500000000 0.2500000000 0.7500000000
S 0.7500000000 0.7500000000 0.2500000000
K_POINTS automatic
12 12 12 0 0 0
CELL_PARAMETERS angstrom
5.4145000000 0.0000000000 0.0000000000
0.0000000000 5.4145000000 0.0000000000
0.0000000000 0.0000000000 5.4145000000

pw.nscf.in is used in the next pw.x calculation:

/opt/intel/oneapi/mpi/2021.1.1/bin/mpirun -np 16 /opt/qe-7.1/bin/pw.x -i

pw.nscf.in > pw.nscf.out &

output:

the Fermi energy is 7.6204 ev

(compare with: 7.6352 eV, computed in scf)

This run was terminated on: 11: 2:52 9Jul2024

=------------------------------------------------------------------------------=
JOB DONE.
=------------------------------------------------------------------------------=

Notice that the Fermi level is slightly lower in this calculation.

5. A density of states calculation can now be performed with the Kohn Sham wavefunctions
having been calculated. To do this, rst write a dos.in le and use it as an input for a dos.x
calculation:

/opt/qe-7.1/bin/dos.x -i dos.zns.in > dos.zns.out &

This run was terminated on: 11: 8:28 9Jul2024

=------------------------------------------------------------------------------=
JOB DONE.
=------------------------------------------------------------------------------=

The data can be plotted with gnuplot, but this requires another le (dos.gp) that controls the
output, axes, etc.

fi
fi
fi
 4

# set Fermi energy to correct value

Efermi=7.620
# ... and uncomment the following line
set yzeroaxis lt -1

set grid
set xlabel "Energy (eV)"
set ylabel "DOS"
set style fill solid 0.4
set format y "%.1f"

set title "Total DOS\n( press <Enter> in the terminal to exit ... )"

plot [:][:] \
'zns.dos' u ($1-Efermi):($1-Efermi<0?$2:0) notit w filledcurve y=0 lt 1, \
'zns.dos' u ($1-Efermi):2 notit w l lt 1 lw 2
pause -1

To plot the data, use gnuplot:

gnuplot dos.gp

Total DOS
( press <Enter> in the terminal to exit ... )
70.0

60.0

50.0

40.0
DOS

30.0

20.0

10.0

0.0

-10.0
-30 -25 -20 -15 -10 -5 0 5 10
Energy (eV)

6. Plot the band structure. This is kind of tricky and will require you to identify lattice points where
the energies will be computed. In brief, unlike an isolated molecule whose orbital energies are
xed, in the solid state, the energies of bands varies across the lattice. For this reason, is is
customary to plot energy bands in terms of particular lattice points. To my knowledge, there is
no algorithm for determing the best path of these points. Rather, you should look at what the
fi

 5

standard is for the crystal class for your material. In this case, I chose 3 lattice points, L, Γ, and
X as those were reported in the original paper on which this demonstration was based.

A. Write a bands.in le where the line following "K points automatic" is substitued for teh
coordinates, in order, of the high-symmetry points in the band path. One way to determine
these coordinates is to open the original .cif in seek path, another tool available from
materialscloud.org.4 As with the pw input generator, use of this tool should be cited
appropriately.5 Among other pieces of information will be a "suggested path" as well as
coordinates (in terms of the lattice structure) of these high symmetry points. I'm ignoring
this suggested path as I have a publication in hand with the path against which these
results will be compared:

Suggested path
Γ—X—U|K—Γ—L—W—X|Γ—X'—U'|K'—Γ—L'—W'—X'

High-symmetry points (scaled units)

Label k1 k2 k3
Γ 0.0000000000 0.0000000000 0.0000000000
K 0.3750000000 0.3750000000 0.7500000000
K' -0.3750000000 -0.3750000000 -0.7500000000
L 0.5000000000 0.5000000000 0.5000000000
L' -0.5000000000 -0.5000000000 -0.5000000000
U 0.6250000000 0.2500000000 0.6250000000
U' -0.6250000000 -0.2500000000 -0.6250000000
W 0.5000000000 0.2500000000 0.7500000000
W' -0.5000000000 -0.2500000000 -0.7500000000
W2 0.7500000000 0.2500000000 0.5000000000

4 https://www.materialscloud.org/work/tools/seekpath
5 If you use this tool, please cite the following work:
• Y. Hinuma, G. Pizzi, Y. Kumagai, F. Oba, I. Tanaka, Band structure diagram
paths based on crystallography, Comp. Mat. Sci. 128, 140 (2017). DOI:
10.1016/j.commatsci.2016.10.015 (the "HPKOT" paper; arXiv version:
arXiv:1602.06402).
• You should also cite Spglib that is an essential library used in the
implementation: A. Togo, I. Tanaka, "Spglib: a software library for crystal
symmetry search", arXiv:1808.01590 (2018)
• The input parsers use a number of libraries (see name in the dropdown list)
from ASE, qe-tools or pymatgen.

fi

 6

W 2' -0.7500000000 -0.2500000000 -0.5000000000

X 0.5000000000 0.0000000000 0.5000000000
X' -0.5000000000 0.0000000000 -0.5000000000

As we are concerned only with L, Γ, and X, these are the only points we'll need in pw.scf.bands.in.
The order is important and the #'s following the coordinates is the number of k points to be
calculated in band structure for each high-symmetry lattice point. As shown, the path goes from L
- Γ - X - L.

&CONTROL
calculation = 'bands'
etot_conv_thr = 8.0000000000d-05
forc_conv_thr = 1.0000000000d-04
outdir = './out/'
prefix = 'aiida'
pseudo_dir = './pseudo/'
tprnfor = .true.
tstress = .true.
verbosity = 'high'
/
&SYSTEM
degauss = 1.4699723600d-02
ecutrho = 7.2000000000d+02
ecutwfc = 9.0000000000d+01
ibrav = 0
nat = 8
nosym = .false.
ntyp = 2
occupations = 'smearing'
smearing = 'cold'
/
&ELECTRONS
conv_thr = 1.6000000000d-09
electron_maxstep = 80
mixing_beta = 4.0000000000d-01
/
ATOMIC_SPECIES
S 32.065 s_pbesol_v1.4.uspp.F.UPF
Zn 65.38 zn_pbesol_v1.uspp.F.UPF
ATOMIC_POSITIONS crystal
Zn 0.0000000000 0.0000000000 0.0000000000
Zn 0.0000000000 0.5000000000 0.5000000000
Zn 0.5000000000 0.0000000000 0.5000000000
Zn 0.5000000000 0.5000000000 0.0000000000
S 0.2500000000 0.2500000000 0.2500000000
S 0.2500000000 0.7500000000 0.7500000000
S 0.7500000000 0.2500000000 0.7500000000
S 0.7500000000 0.7500000000 0.2500000000
K_POINTS {crystal_b}
4



 7

# L-Gamma-X-L
0.5000000000 0.5000000000 0.5000000000 20
0.0000000000 0.0000000000 0.0000000000 20
0.5000000000 0.0000000000 0.5000000000 20
0.5000000000 0.5000000000 0.5000000000 0
CELL_PARAMETERS angstrom
5.4145000000 0.0000000000 0.0000000000
0.0000000000 5.4145000000 0.0000000000
0.0000000000 0.0000000000 5.4145000000

B. Calculate the eigenvalues at the k-points along the speci ed path:

/opt/qe-7.1/bin/pw.x < pw.zns.bands.in > pw.zns.bands.out &

C. Write a zns.bands.in le that takes the output above and renders a .dat le for plotting.
Apply bands.x to this input le:

&BANDS
prefix = 'aiida'
outdir = './out/'
filband = 'zns.bands.dat'
lsym = .true.,
/

/opt/qe-7.1/bin/bands.x < bands.zns.in > bands.zns.out

D. Write a le (spaghetti.gp), that will read the bands.zns.dat le and generate a band plot.
Insert the correct letter/symbol for high-symmetry points used in your path. The Fermi level
is obtained from pw.nscf.out. The x-coordinate of these values is obtained from the
bands.zns.out le as shown:

Reading collected, re-writing distributed wavefunctions

high-symmetry point: 0.5000 0.5000 0.5000 x coordinate 0.0000
high-symmetry point: 0.0000 0.0000 0.0000 x coordinate 0.8660
high-symmetry point: 0.5000 0.0000 0.5000 x coordinate 1.5731
high-symmetry point: 0.5000 0.5000 0.5000 x coordinate 2.0731

Plottable bands (eV) written to file zns.bands.dat.gnu

Bands written to file zns.bands.dat

Band structure of ZnS without Hubbard Correction


fi
fi
fi
fi
fi
fi
fi
 8
8.0

6.0

4.0

2.0
Energy (Ry)

0.0

-2.0

-4.0

-6.0

-8.0
L Γ X L

spaghetti.gp

# set Fermi energy to correct value

Efermi=7.620
# ... and uncomment the following line
set xzeroaxis lt -1

set grid xtics lt -1 lw 1

set format y "%5.1f"
set format x ""
set ylabel "Energy (Ry)"
unset xlabel

set xtics ("L" 0.0000, "{/Symbol G}" 0.8660, "X" 1.5731, "L" 2.0731)

plot [0:2.0731] 'zns.bands.dat.gnu' u 1:($2-Efermi) notitle w linespoints

lw 2 p
t 0
pause -1

E. Plot the band structure:

gnuplot spaghetti.gp

The bandgap, which is direct, is the 1.95 eV separation at the Γ point. The Fermi level is set at 0 in
this analysis.


 9

Hubbard Correction

But, But, But, DFT is known to underestimate the bandgap in materials and to the point where
semiconductors can appear metallic.6 To correct for this, the DFT + U correction was introduced
by Anisimov7 and developed by others. To do this, a U parameter is selected by experience or
optimization of the wavefunction with systematic variation in U. In this example, we consider the
former applied to zincblende.

Method

Modify .nscf.in to include hubbard parameters (here, they are included after the K_POINTS
section).8

HUBBARD {ortho-atomic}
U Zn-3d 13.7
U S-3p 3.51

Next, run the calculation through including to the band calculation including the HUBBARD
parameters to obtain a new band structure (shown below). The new bandgap is 3.05 eV. Additional
corrections further improve the match with experiment.

Band stucture of ZnS with Hubbard correction

6 Dudarev, S. L., Botton, G. A., Savrasov, S. Y., Humphreys, C. J. & Sutton, A. P. Electron-energy-
loss spectra and the structural stability of nickel oxide: An lsda+u study. Phys. Rev. B 1998 57,
1505–1509.
7Anisimov, V. I., Zaanen, J. & Andersen, O. K. Band theory and Mott insulators: Hubbard U instead
of Stoner I. Phys. Rev. B 44, 943–954 (1991).
8 Selected from reference 1.

 10

Additional citations:

ASE: Ask Hjorth Larsen, Jens Jørgen Mortensen, Jakob Blomqvist, Ivano E. Castelli, Rune
Christensen, Marcin Dułak, Jesper Friis, Michael N. Groves, Bjørk Hammer, Cory Hargus, Eric D.
Hermes, Paul C. Jennings, Peter Bjerre Jensen, James Kermode, John R. Kitchin, Esben Leonhard
Kolsbjerg, Joseph Kubal, Kristen Kaasbjerg,Steen Lysgaard, Jón Bergmann Maronsson, Tristan
Maxson, Thomas Olsen, Lars Pastewka, Andrew Peterson, Carsten Rostgaard, Jakob Schiøtz, Ole
Schütt, Mikkel Strange, Kristian S. Thygesen, Tejs Vegge, Lasse Vilhelmsen, Michael Walter,
Zhenhua Zeng, Karsten Wedel Jacobsen The Atomic Simulation Environment—A Python library for
working with atoms
J. Phys.: Condens. Matter 2017, 29, 273002.

S. R. Bahn and K. W. Jacobsen An object-oriented scripting interface to a legacy electronic

structure code Comput. Sci. Eng., 2002, 4, 56-66, 2002

qe-tools: S. P. Huber et al., AiiDA 1.0, a scalable computational infrastructure for automated
reproducible work ows and data provenance, Scienti c Data 2020, 7, 300. DOI: 10.1038/
s41597-020-00638-4

M. Uhrin et al., Work ows in AiiDA: Engineering a high-throughput, event-based engine for robust
and modular computational work ows, Computational Materials Science 2021, 187, 110086; DOI:
10.1016/j.commatsci.2020.110086

pymatgen: Shyue Ping Ong, William Davidson Richards, Anubhav Jain, Geo roy Hautier,
Michael Kocher, Shreyas Cholia, Dan Gunter, Vincent Chevrier, Kristin A.
Persson, Gerbrand Ceder. *Python Materials Genomics (pymatgen) : A Robust,
Open-Source Python Library for Materials Analysis.* Computational Materials
Science, 2013, 68, 314–319. https://doi.org/10.1016/j.commatsci.2012.10.028

fl
fl
fl
fi
ff