Structural and electronic properties of semiconductors and metals
 Prev: LabQSM#Module 1: Basic DFT calculations and Convergences (9h)
 Next: Non selfconsistent calculations: Band structures and Density Of States
Structural and electronic properties of Diamond
In this tutorial we will see how to setup a calculation and to get total energies using the PW code from the Quantum ESPRESSO distribution.
Some helpful conversions:
1 Bohr = 1 a.u. (atomic unit) = 0.529177249 Angstroms. 1 Rydberg = 13.6056981 eV 1 eV =1.60217733 x 10^19 Joules
For all firstprinciples calculations, you must pay attention to two convergence parameters. The first one is the energy cutoff, which is the max kinetic energy used in wavefunction expansion. The second is the number of kpoints, which measures how well the continuous integral over the BZ is discretized.
Diamond is a facecentered cubic structure with two C atoms at 0 0 0 and 0.25 0.25 0.25 a is the lattice parameter
Now let's see in detail how a QE input is structured to make a total energy calculation for this system. An example can be found in ~/LabQSM/LAB_1/test_diamond/scf.diamond.in
that you can read e.g. using the editor vi
:
Input file description
&CONTROL prefix='diamond', calculation = 'scf' restart_mode='from_scratch', pseudo_dir = './pseudo/' outdir = './SCRATCH' / &SYSTEM ibrav = 2, celldm(1) = 6.7402778 nat = 2, ntyp = 1, ecutwfc = 40.0 / &ELECTRONS mixing_mode = 'plain' mixing_beta = 0.7 conv_thr = 1.0d8 / ATOMIC_SPECIES C 12.011 C.pzvbc.UPF ATOMIC_POSITIONS (alat) C 0.00 0.00 0.00 C 0.25 0.25 0.25 K_POINTS {automatic} 4 4 4 0 0 0
The input file for PWscf is structured in a number of NAMELISTS and INPUT CARDS:
&NAMELIST1 ... / &NAMELIST2 ... / &NAMELIST3 ... / INPUT_CARD1 .... .... INPUT_CARD2 .... ....
The use of NAMELISTS allows to specify the value of an input variable, not all the variables need to be specified, for most variable, a default value is assigned. Variable can be inserted in any order. NAMELISTS are read in a specific order. INPUT CARDS are specific of QuantumESPRESSO codes and are used to provide input data that are always needed.
There are three mandatory NAMELISTS (while more can be required according to the calculation type):
 &CONTROL: variables that control the kind of calculation (here scf), where to get the pseudopopotentials, and the verbosity needed in the output.
 &SYSTEM: variables specifying the system as the crystal structure, number of atoms, dimension of the basis set.
 &ELECTRONS: variables controlling the algorithm to solve the KohnSham equations.
 &IONS: needed only for some values of the
calculation
variable (CONTROL namelist), such as set to "relax" or "vcrelax".  &CELL: as above, needed for selected kind of calculations, including "vcrelax".
Next, we have three mandatory INPUT CARDS:
 ATOMIC_SPECIES : name, mass and pseudopotential file
 ATOMIC_POSITIONS: atom type and coordinate in the unit cell
 K_POINTS: definition of the k point grid for the BZ integration
Other NAMELISTS relative to other kinds of calculations as relaxations, cell relaxations and other INPUT CARDS will be shown later.
The complete list of QE variables can be found in the documentation of the code at:  http://www.quantumespresso.org/Doc/INPUT_PW.html
Now that we know how a QE input looks like let's run it and see what the code does and analyze the output:
Running pw.x
Let's create a working directory e.g.
$> mkdir diamond
and copy in this directory the scf.diamond.in input file we have inspected.
Now you can run the pw.x executable as:
$> pw.x < scf.diamond.in > scf.diamond.out
It is also possible to run in parallel:
$> export OMP_NUM_THREADS=2 #we use 2 OpenMP threads $> pw.x < scf.diamond.in > scf.diamond.out
or
$> mpirun np 2 pw.x < scf.diamond.in > scf.diamond.out #we use 2 MPI tasks running on 2 processors
Several files will be generated in ./SCRATCH
, most of them needed mainly for postprocessing now let's have a look to the readable output (scf.diamond.out
).
Output file description
Bearing in mind the logical flow of a selfconsistent KohnSham DFT calculation performed using the density mixing approach, defined in the figure,
we can now inspect the output (scf.diamond.in > scf.diamond.out) using the vi/more/less
command.
This is what we find:
A recap of the configuration that is being calculated:
bravaislattice index = 2 lattice parameter (alat) = 6.7403 a.u. unitcell volume = 76.5550 (a.u.)^3 number of atoms/cell = 2 number of atomic types = 1 number of electrons = 8.00 number of KohnSham states= 4 kineticenergy cutoff = 40.0000 Ry charge density cutoff = 160.0000 Ry convergence threshold = 1.0E08 mixing beta = 0.7000 Exchangecorrelation= SLA PZ NOGX NOGC
Detailed information of the unit cell:
celldm(1)= 6.740278 celldm(2)= 0.000000 celldm(3)= 0.000000 celldm(4)= 0.000000 celldm(5)= 0.000000 celldm(6)= 0.000000
crystal axes: (cart. coord. in units of alat) a(1) = ( 0.500000 0.000000 0.500000 ) a(2) = ( 0.000000 0.500000 0.500000 ) a(3) = ( 0.500000 0.500000 0.000000 )
reciprocal axes: (cart. coord. in units 2 pi/alat) b(1) = ( 1.000000 1.000000 1.000000 ) b(2) = ( 1.000000 1.000000 1.000000 ) b(3) = ( 1.000000 1.000000 1.000000 )
Information on the used pseudopotential:
PseudoPot. # 1 for C read from file: /home/max/LabQSM/pseudo/C.pzvbc.UPF MD5 check sum: ab53dd623bfeb79c5a7b057bc96eae20 Pseudo is Normconserving, Zval = 4.0 Generated by new atomic code, or converted to UPF format Using radial grid of 269 points, 1 beta functions with: l(1) = 0
atomic species valence mass pseudopotential C 4.00 12.01100 C ( 1.00)
Information on the symmetries operation detected
48 Sym. Ops., with inversion, found (24 have fractional translation)
Information on the k points used to sample the Brillouin zone and FFT gird:
number of k points= 8 cart. coord. in units 2pi/alat k( 1) = ( 0.0000000 0.0000000 0.0000000), wk = 0.0312500 k( 2) = ( 0.2500000 0.2500000 0.2500000), wk = 0.2500000 k( 3) = ( 0.5000000 0.5000000 0.5000000), wk = 0.1250000 k( 4) = ( 0.0000000 0.5000000 0.0000000), wk = 0.1875000 k( 5) = ( 0.7500000 0.2500000 0.7500000), wk = 0.7500000 k( 6) = ( 0.5000000 0.0000000 0.5000000), wk = 0.3750000 k( 7) = ( 0.0000000 1.0000000 0.0000000), wk = 0.0937500 k( 8) = ( 0.5000000 1.0000000 0.0000000), wk = 0.1875000
Dense grid: 2685 Gvectors FFT dimensions: ( 20, 20, 20)
Here note that the 2685 Gvectors are defined, using the ecutrho variable, as points inside a cutoff sphere in reciprocal space. The FFT grid, mapping direct to reciproca space, is built as a box containing the cutoff sphere.
Intermediate energies calculated during the scf loop
Selfconsistent Calculation
iteration # 1 ecut= 40.00 Ry beta= 0.70 Davidson diagonalization with overlap ethr = 1.00E02, avg # of iterations = 2.0
total cpu time spent up to now is 0.2 secs
total energy = 22.67865569 Ry HarrisFoulkes estimate = 22.74103732 Ry estimated scf accuracy < 0.12254037 Ry ...
Near the end “hopefully” our converged results
Final Eigenvalues for each k point:
End of selfconsistent calculation
k = 0.0000 0.0000 0.0000 ( 331 PWs) bands (ev):
8.0812 13.3868 13.3868 13.3868
k =0.2500 0.25000.2500 ( 323 PWs) bands (ev):
6.3642 6.7040 11.6380 11.6380
For each kpoint, the number of Gvectors used to represent the wave functions are provided
(331, 323, ...). These are kdependent (the cutoff sphere for wave function is defined on k+G vectors and is therefore kdependent).
When ecutrho = 4* ecutwfc
(as it should for notaconserving pseudopotentials), the number of Gvectors in the density grid is almost 8 times
that of Gvectors in the wfc grid (geometrically, the kinetic energy cutoff is the square of the sphere radius, the number of Gvectors being proportional to its volume).
Final total energy and contribution of the various terms:
Note that the converged total energy is preceded by an exclamation mark "!".
! total energy = 22.68935940 Ry HarrisFoulkes estimate = 22.68935940 Ry estimated scf accuracy < 1.5E09 Ry
The total energy is the sum of the following terms:
oneelectron contribution = 8.08269762 Ry hartree contribution = 1.85703398 Ry xc contribution = 7.05484267 Ry ewald contribution = 25.57424832 Ry
convergence has been achieved in 6 iterations
At the end of the file, you can find info on timings spent in each part of the code (subroutines):
PWSCF : 0.13s CPU 0.32s WALL
Now we can, for instance, extract the total energy calculated at each step of the selfconsistent loop:
$>grep e "total energy" scf.diamond.out
or the accuracy reached AT each step:
$> grep e "scf" scf.diamond.out
or we can grep both of them and put in a twocolumn format:
$> grep e "total energy" e "scf" scf.diamond.out  awk '/l e/{e=$(NF1)}/ scf / {print e, $(NF1)}' 22.67865569 0.12254037 22.68899239 0.00212159 22.68934113 0.00007686 22.68935832 0.00000263 22.68935938 0.00000004 22.68935940 1.5E09
or we can redirect to a file to be plot
$> grep e "total energy" e "scf" scf.diamond.out  awk '/l e/{e=$(NF1)}/ scf / {print e, $(NF1)}' > diamond_scf.dat $> gnuplot $> plot "diamond_scf.dat" u 0:1 wlp $> p "diamond_scf.dat" u 0:(log10($2)) w lp ,8

Total Energy vs iteration

Relative error vs iteration
Convergences
Exercise 1
Parsing pw output files
 Write a script to extract the total energy and the lattice parameter from a pw output file
 Print them on the same output line
 Solve the problem using a bash script (e.g. called
extract.sh
)  Try to do the same also using a awk script (
extract.awk
)  If you are a python expert you can give it a try (
extract.py
)
Note that .sh and .py scripts can be made easily working on multiple files
 Hint1: the total energy is marked by "!", while the lattice parameter can be taken from celldm1 (marked by "lattice parameter")
 Hint2: look at an example and build the scripts checking step by step, use e.g. scf.diamond.out
Exercise 2
Convergence wrt kpoints
 Using a script, perform different runs in order to converge the absolute energy of your system with respect to the kpoint sampling
(keep other variables fixed).
 You can use a cutoff lower than the converged one (~40 Ry)
 Hint: Try 4x4x4, 6x6x6, 8x8x8 etc...
 Take also a look at the number of the irreducible kpoints generated by the code and at the relative execution time.
Exercise 3
Convergence wrt wavefunction kinetic energy cutoff (ecutwfc)
 Using a script, perform different runs in order to converge the absolute energy of your system with respect to the plane waves kinetic energy cutoff used to represent the wave functions.
 Hint: Use an increment 20 Ry in the range 20200 Ry;
 for the purpose of evaluating the scaling of timetosolution wrt ecutwfc, you can push the calculations to even larger (and further unrealistic) values (such as 300 and 400 Ry).
 Convergence criteria ~5meV/atom
BTW: Note that this is NOT the procedure to follow to determine the kinetic energy cutoff to use in production runs.
In fact, values determined in this way tend to be much larger than the cutoff actually needed for most applications we are interested in.
Exercise 4
Convergence of Energy differences wrt ecutwfc
 Using a script, compute the energy difference corresponding to using two lattice parameters and perform several runs in order to converge such a difference with respect to the plane wave kinetic energy cutoff.
 Hint1: Use an increment of 20 Ry in the range 20200 Ry
 Hint2: consider a variation of the lattice parameter providing a sensible energy difference (12% is ok)
 Convergence criteria: ~ 0.1 mRy/atom
 How does the convergence compare with the one of the bare total energy ?
Exercise 5
Convergence of forces wrt ecutwfc
 Displace a C atom by +0.05 (in alat or crystal relative coords) in the z direction to induce nonzero forces acting atoms (otherwise forbidden by symmetry).
 Keeping the other parameters fixed, calculate the forces on C as a function of the kinetic energy cutoff.
 A good threshold on force convergence is 1 to 0.1 mRy/Bohr.
Hint: Remember to set tprnfor=.true.
in the namelist control.
Lattice parameter & Elastic constants
Here we briefly introduce some poorman approaches to compute the lattice parameter and bulk modulus of simple crystals by directly evaluating total first and second derivatives wrt the lattice parameter.
We do this in the form of guided exercises.
Exercise 6
Lattice parameter of diamond
 Using a script, perform different runs in order to obtain the lattice parameter of Diamond.
 Keep fixed the kinetic energy cutoff and kpoint sampling to the values converged before for energy differences (or run the calculations for increasing values of ecutwfc).
 Hint: strain the experimental lattice parameter from 3% to +3% with steps of 1%
 Exp: Alat_C = 6.741 Bohr
Exercise 7
Computing the bulk modulus
 Using a script, perform runs at different values of cell volume (ie using different lattice parameters) and compute the bulk modulus
 This is defined according to the following expressions:
P =  dE/dV B =  Vo dP/dV = Vo d2E/dV2 Vo: equilibrium volume
 Hint1: strain the experimental lattice parameter from 3% to +3% with steps of 1%, collect results and fit them against volume.
 Perform the second derivative numerically by fitting E=E(V) (see folder
LabQSM/tools
for info about how to do this nuerically)  The python script
LabQSM/tools/compute_B_modulus.py
implements the procedure and provides the constants for unit conversion from a.u. to SI.  Hint2: The following can help with units: [1]
 Note1: primitive cell volume = 1/4 conventional cell ;
 Note2: Exp value for diamond: 443 GPa; See also R. Gaudoin and W.M.C. Foulkes, Phys. Rev. B 66, 052104 (2002).
Structural Relaxation
This is a new runlevel of Quantum ESPRESSO and it is aimed at mechanical stability, i.e. making total energy stationary wrt FORCES and components of the STRESS TENSOR.
 This is a necessary and sufficient condition for equilibrium of a periodic array of atoms; note that zero FORCES alone are not sufficient.
 In practice, the code will perform a series of scf calculation moving ions or cell parameters until a convergence criterium is reached.
These kinds of calculation are activated selecting:
calculation=”relax” # cell is fixed
Or
calculation=”vcrelax” # variablecell relax
and they need additional mandatory namelists which are:
 &IONS # cell is fixed, but atoms move: it will contain input variables that control the ionic motion dynamic in structural relaxation
 &CELL # the cell moves: it will contain variables that control the shellshape dynamics in structural relaxation
Examples are:
&CONTROL [...] forc_conv_thr=1.0d4 etot_conv_thr=1.0d5 / &SYSTEM [...] / &ELECTRONS [...] / &IONS ion_dynamics=”bfgs” [...] / &CELL cell_dynamics=”bfgs” press_conv_thr=0.5 cell_dofree=”all” [...] /
In this example, we used the BFGS quasinewton algorithm for both cell dynamics and ion dynamics, and we are optimizing all generating vectors and angles (cell_dofree=”all”).
It is also possible to keep fixed some structural parameter fixed or insert constraints see cell_dofree options.
Exercise 8
Diamond relaxation:
 Take a scf input file for diamond, offset one of the two atoms
 Set a initial lattice far from the equilibrium condition (eg 5%)
 Using “relax” and “vcrelax” find the optimum geometry
 Plot the resulting lattice parameter as a function of the kinetic energy cutoff
Hint: Pay attention to the units of cell dimension and atomic position (as they are the in initial alat unit)
Note that vcrelax is particularly sensitive to the kinetic energy cutoff
A metallic system: Aluminium
Aluminium is even simpler than Diamond as we have just one atom per unit cell in a fcc lattice. Important: it is a metal, we cannot consider valence bands only and we need to sample the Fermi surface. To do that we need to "smear" the occupation numbers (e.g. with a gaussian smearing). In the &SYSTEM namelist we include the following variables.
occupations=’smearing’ smearing=’gaussian’ degauss=0.01
Othe forms of the smearing function are available, and you can found them in the documentation
Exercise 9
Convergences in a metallic system:
 build an scf input to calculate total energy for Aluminium (pseudopotential can be found in the usual /home/max/LabQSM/pseudo directory), experimental lattice parameter: 4.0495 \AA
 Build a script to check the convergence with respect k point sampling and smearing (note that the two variables are not independent)
 Plot the total energy with respect the smearing value for the different set of samplings
Hint: use values nk=4,8,12,16 and a set of smearing in the interval 0.10.01 Ry