Self consistent field calculation for silicon
We need to provide various important parameters for the self consistent
calculation (solves the KohnSham equation selfconsistently) via an input file.
In QE input files, there are NAMELISTS
and INPUT_CARDS
. NAMELISTS
variables have default values, and new values can be provided as required for a
specific calculation. The variables can be declared in any specific order. On
the other hand, the variables in the INPUT_CARDS
has always to be specified
and in specific order. Logically independent INPUT_CARDS
may be organized in
any order.
There are three mandatory NAMELISTS
in PWscf
: (1) &CONTROL
: specifies the
flux of computation, (2) &SYSTEM
: specifies the system, and (3) &ELECTRONS
:
specifies the algorithms used to solve the KohnSham equation. There are two
other NAMELISTS
: &IONS
and &CELLS
, which need to be specified depending on
the calculation.
Three INPUT_CARDS
: ATOMIC_SPECIES
, ATOMIC_POSITIONS
, and K_POINTS
in
PWscf
are mandatory. There are few others that must be provided in certain
calculations.
Below is our input file pw.scf.silicon.in for silicon in standard diamond (FCC)
structure. The input files are typically named with .in
prefix, while output
files are named with .out
prefix for their easier identification. The input
parameters are organized in &namelists
followed by their fields or cards. The
&control
, &system
, and &electrons
namelists are required. There are also
optional &cell
and &ions
, you must provide them if your calculation require
them. Most parameters in the namelists
have default values (which may or may
not suit your needs), however some variables you must always provide. Comment
lines can be added with lines starting with a !
like in Fortran. Also,
parameter names are not casesensitive, i.e., &control
and &CONTROL
are the
same.
&CONTROL
! we want to perform self consistent field calculation
calculation = 'scf',
! prefix is reference to the output files
prefix = 'silicon',
! output directory. Note that it is deprecated.
outdir = './tmp/'
! directory for the pseudo potential directory
pseudo_dir = '../pseudos/'
! verbosity high will give more details on the output file
verbosity = 'high'
/
&SYSTEM
! Bravais lattice index, which is 2 for FCC structure
ibrav = 2,
! Lattice constant in BOHR
celldm(1) = 10.26,
! number of atoms in an unit cell
nat = 2,
! number of different types of atom in the cell
ntyp = 1,
! kinetic energy cutoff for wavefunction in Ry
ecutwfc = 30
! number of bands to calculate
nbnd = 8
/
&ELECTRONS
! Mixing factor used in the selfconsistent method
mixing_beta = 0.6
/
ATOMIC_SPECIES
Si 28.086 Si.pzvbc.UPF
ATOMIC_POSITIONS (alat)
Si 0.0 0.0 0.0
Si 0.25 0.25 0.25
K_POINTS (automatic)
6 6 6 0 0 0
I am using the pseudo potential file (Si.pzvbc.UPF
) downloaded from Quantum
Espresso Website.
You must read the PWscf user manual for indepth understanding. Check the
qex.x/PW/Doc/
folder under your installation directory. There is also another
file INPUT_PW.html
regarding the details of input parameters. PW stands for
plane waves.
Run pw.x
in self consistent mode for silicon.
pw.x < pw.scf.silicon.in > pw.scf.silicon.out
# For parallel execution
mpirun np 4 pw.x inp pw.scf.silicon.in > pw.scf.silicon.out
Note that I have added the executable path to my bash/zsh profile, otherwise you
have to provide the full path where the pw.x
executable is located.
Now let’s look at the output file pw.scf.silicon.out
and see how the
convergence is reached:
grep e 'total energy' e estimate pw.scf.silicon.out
and you should see something like this:
total energy = 15.85014573 Ry
HarrisFoulkes estimate = 15.86899637 Ry
estimated scf accuracy < 0.06093037 Ry
total energy = 15.85194177 Ry
HarrisFoulkes estimate = 15.85292281 Ry
estimated scf accuracy < 0.00462014 Ry
total energy = 15.85218359 Ry
HarrisFoulkes estimate = 15.85220235 Ry
estimated scf accuracy < 0.00011293 Ry
! total energy = 15.85219789 Ry
HarrisFoulkes estimate = 15.85219831 Ry
estimated scf accuracy < 0.00000099 Ry
The total energy is the sum of the following terms:
It is important to note that the absolute value of DFT total energy is not with respect to the vacuum reference, and depends on the chosen pseudopotential. The meaningful measure is the difference in total energy, where various offsets cancel out.
In the above calculation, if you check the output file pw.scf.silicon.out
, you
will find: highest occupied, lowest unoccupied level (eV): 6.2117 6.8442.
Therefore, the bandgap is 0.6325 eV, which is an underestimation of actual
bandgap (1.12 eV).

Reduce
mixing_beta
value, especially if there is an oscillation around the convergence energy. 
If it is a metallic system, use smearing and degauss. In this case, the SCF accuracy gradually goes down then suddenly increases (due to slight change in Fermi energy highest occupied/ lowest unoccupied levels change).

Increase energy and charge density cutoffs (make sure they are sufficient).

Certain pseudo potential files have issues, you may try with pseudo potentials from different libraries.

Suggested values for the
conv_thr
: for energy and eigenvalues (scf calculation) 1.0d7, for forces (relax calculation) 1.0d8, for stress (vcrelax calculation) 1.0d9 Ry. For certain calculation convergence might be very slow for the first iteration, one can start the calculation with a higher threshold, after few iterations reduce it and restart the calculation.
There are several other important information is printed on the output file. Exchange correlation used in the calculation:
Exchangecorrelation= SLA PZ NOGX NOGC
Where SLA
→ Slater exchange; PZ
→ PerdewZunger parametrization of the LDA;
NOGX
and NOGC
indicates that density gradients are not taken into account.
We can see the total number of plane waves (1067) uses in our calculation:
Parallelization info

sticks: dense smooth PW Gvecs: dense smooth PW
Min 108 108 34 1489 1489 266
Max 109 109 35 1492 1492 267
Sum 433 433 139 5961 5961 1067
Number of KohnSham states:
number of electrons = 8.00
number of KohnSham states= 8
In our calculation we have specified the number of bands = 8. Otherwise, there would be 4 bands for 8 electrons in case of non spinpolarized systems.
Resources
 https://www.quantumespresso.org/Doc/pw_user_guide/
 Quantum Espresso Input Generator (can help crating QE input files)