sc :: elementary

# Electronic structure of solids with NwChem

### Introduction

For this exercise we will try to reproduce a band structure and density of states of TiO2 in rutile lattice structure. We will rely on the previously published results of Landmann et al. [1]. We are going to use NwChem program to do the calculation.

Titanium oxide (TiO2) exists in three lattice structures, i.e. rutile, ilmenite and anatase. We will consider the rutile structure that is shown on the picture below. The picture has been prepared using program VESTA [4]. Blue atoms represent titanium and red atoms represent oxygen.

NwChem's module for solid states is called NWPW and is a pseudopotential plane-wave density functional theory code with three main modules i.e. pspw, band and paw. All three modules have different capabilities. Band structure and density of states are possible to do using band program.

#### Bibliography

[1] M.Landmann, E.Rauls and W.G.Schmidt, J.Phys.C.M. 24, 195503 (2012), DOI

[2] W.Setyewan and S.Cutrarolo, Comp. Mat. Sci. 48, 299 (2010), DOI

[3] Pseudopotential plane-wave density functional theory (NWPW) NwChem manual page

[4] K.Momma and F.Izumi, J. Appl. Crystallogr., 44, 1272 (2011) DOI

### Band structure

In NwChem band structure calculation can be done using a subprogram band. The same program calculates the density of states. band can run several tasks in NwChem e.g. optimize, energy, structure, dos and band_dplot. Listed tasks are responsible for distinct types of calculations i.e. geometry optimization, energy calculation (single point calculation), band structure calculation, density of states, generation data for visualization (e.g. density, orbitals) respectively.

All keywords for band program have to be present in the nwpw block

nwpw
...
end


where jobType is one of optimize, energy, structure, dos or band_dplot.

First step in the calculation of band structure is the correct specification of the crystal unit cell and the coordinates of atoms in the unit cell.

#### Geometry of the crystal

In our example we will examine titanium oxide TiO2 in the rutile structure. We will specify the dimensions of the unit cell in Angstroms and angles in degrees. TiO2 in rutile structure has tetragonal unit cell $\alpha=\beta=\gamma=90^\circ$ and lattice parameters $a=b \neq c$. The symmetry of the rutile structure is represented by $P42/mnm$ group. The exact lattice structure parameter for our example will be taken from the online supplement to "Structure of materials" book TiO2 data as a=b=4.5937A and c=2.9581A. The coordinates of atoms are given in the fractional units of the lattice structure parameters. The full geometry block looks like this

geometry units angstroms center noautosym noautoz print
system crystal
lat_a 4.5937
lat_b 4.5937
lat_c 2.9581
alpha 90.0
beta  90.0
gamma 90.0
end
symmetry P42/mnm
Ti  0.0000 0.0000 0.0000
O   0.3053 0.3053 0.0000
end


We have specified only symmetry unique atoms in the unit cell and explicitly specified the symmetry using symmetry P42/mnm statement. The same can be achieved by explicitly giving coordinated of 6 atoms in the unit cell and skipping the symmetry keyword.

#### Density Functional Theory for the energy

The calculation starts with the calculation of the DFT energy. The energy task input may read

nwpw
ewald_rcut 3.0
ewald_ncut 8
lmbfgs
xc pbe96
np_dimensions -1 -1 4
monkhorst-pack 9 9 9
end


In the above example we will use PBE96 functional xc pbe96 and lmbfgs solver for the energy calculation.

Keyword Description
ewald_rcut cutoff radius in Ewald summation
ewald_ncut number of unit cells in the real space Ewald summation
lmbfgs Limited-memory BFGS solver
xc pbe96 functional from the list of Vosko (default), LDA, PBE96, revPBE, PBEsol, HSE
np_dimensions FFT, Orbital, K-space (-1 sets value to default)
monkhorst-pack Monkhorst-Pack sampling of the Brillouin zone

#### Brillouin Zone scan

band program has an automatic scan over the high symmetry points implemented as a walk through the Brillouin zone via set of specific K-points.

The definition of the Brillouin zone scan in NwChem input file is a block starting with keyword brillouin_zone and ending with end which is encapsulated in the nwpw block. The general input will look like this

brillouin_zone
zone_name UserDefinedName
path LatticeType K-Points-list
end


allowed keywords for LatticeType and the list of allowed K-Points-list for that lattice type are listed in the table below.

Structure Keyword K-points
Simple cubic - SC sc gamma, m, r, x
Face centered cubic - FCC fcc gamma, k, l, u, w, x
Body centered cubic - BCC bcc gamma, h, n, p
Rhombohedral rhombohedral not currently implemented as of NwChem 6.5
Hexagonal hexagonal gamma, a, h, k, l, m
Simple tetragonal simple-tetragonal gamma, a, m, r, x, z
Simple orthorhombic simple-orthorhombic gamma, r, s, t, u, x, y, z
Body centered tetragonal bct gamma, m, n, p, x

In the example of TiO2 in rutile simple tetragonal structure we will use the path through the K-points in order suggested by Setyewan and Cutrarolo [2] i.e. path simple-tetragonal gamma x m gamma z r a. The full block specifying the Brillouin zone path reads

brillouin_zone
zone_name tetragonalzone
path simple-tetragonal gamma x m gamma z r a
end


One is not restricted to using the path keyword and a manual specification of the K-points is also possible. In this introductory example we will use the automatic walk via path keyword. Keyword zone_name in the above example is a user defined name for the Brillouin zone calculation. The Brillouin zone scan requires specifying how many virtual orbitals will be included in the calculation. This will define how many bands are calculated.

#### Complete input file

The complete NwChem input file for the band structure calculation using band program may look like this.

echo
title "Rutile band structure calculation - example"

memory 5000 mb

permanent_dir ./
scratch_dir ./

start rutile_example

geometry units angstroms center noautosym noautoz print
system crystal
lat_a 4.5937
lat_b 4.5937
lat_c 2.9581
alpha 90.0
beta  90.0
gamma 90.0
end
symmetry P42/mnm
Ti  0.0000 0.0000 0.0000
O   0.3053 0.3053 0.0000
end

nwpw
ewald_rcut 3.0
ewald_ncut 8
lmbfgs
xc pbe96
np_dimensions -1 -1 4
monkhorst-pack 9 9 9
end

nwpw
virtual 40
brillouin_zone
zone_name tetragonalzone
path simple-tetragonal gamma x m gamma z r a
end
zone_structure_name tetragonalzone
end


The input has four main sections 1) hardware specific keywords like memory,permanent_dir and scratch_dir, 2) crystal structure definition (discussed above), 3) energy calculation (discussed above) and 4) Brillouin zone path scan for the band structure.

In the input file above we gave the name rutile_example to the job and all of the important files will contain that name together with an extension. After the job has completed you should have file with name rutile_example.restricted_band in the permanent folder specified at the top of the input file.

### Density of states

The band program can only evaluate the total density of states. In our example we will run it as a separate job from the Brillouin zone scan. The last part of the input file will differ and will contain the taks band dos input.

set dos:npoints 1000
nwpw
virtual 40
dos-grid 11 11 11
end


The first line set dos:npoints 1000 sets density of the energy grid for density of state. Two following lines i.e. virtual 40 and dos-grid 11 11 11 set the number of virtual orbitals included in the calculation (setting the upper energy limit for the DOS analysis) and the grid size for DOS calculation.

The complete input file may look like this

echo
title "Rutile density of states calculation - example"

memory 5000 mb

permanent_dir ./
scratch_dir ./

start rutile_example

geometry units angstroms center noautosym noautoz print
system crystal
lat_a 4.5937
lat_b 4.5937
lat_c 2.9581
alpha 90.0
beta  90.0
gamma 90.0
end
symmetry P42/mnm
Ti  0.0000 0.0000 0.0000
O   0.3053 0.3053 0.0000
end

nwpw
ewald_rcut 3.0
ewald_ncut 8
lmbfgs
xc pbe96
np_dimensions -1 -1 4
monkhorst-pack 9 9 9
end

set dos:npoints 1000
nwpw
virtual 40
dos-grid 11 11 11
end


The density of states analysis will be saved to separate file i.e. rutile_example.dos.

### Results

The complete picture that we will try to get is present in the outputs of two calculations and two files rutile_example.restricted_band and rutile_example.dos. To have the analysis of the results complete we will also use the output of the first calculation.

#### Labeling Brillouin zone points

The band structure plot can be done directly using the rutile_example.restricted_band file. However, labeling of the x-axis will require analysis of the K-points and matching them with the path length number. The rutile_example.restricted_band contains path length in the first column and energy of bands in the following columns.

To label the x-axis with the characteristic K-points of given lattice we need to search for all high symmetry K-points points on the path in the output file. In our case the selected path was path simple-tetragonal gamma x m gamma z r a. We need to look for the K-points labeled $\Gamma$, $X$, $M$, again $\Gamma$, $Z$, $R$ and $A$. These points are not marked in the output file with their symmetry label, instead, coordinates of the K-point are given. We need to look for the coordinates of the symmetry points in this order $(0,0,0)\,\Gamma$, $(0,0.5,0)\,X$, $(0.5,0.5,0)\,M$, again $(0,0,0)\,\Gamma$, $(0,0,0.5)\,Z$, $(0,0.5,0.5)\,R$, $(0.5,0.5,0.5)\,A$ (consult the picture below on the position of the K-points in the Brillouin zone of tetragonal lattice). Each of these point will have assigned path length that points to its position on the x-scale of the band structure plot.

#### Complete result

On the picture below we present the complete results from this calculation plotted side-by-side. The left panel shows the band structure of TiO2 in rutile structure and the right panel shows the density of states. On both panels the y axis marks the energy in eV.

Energy has been shifted to mark the highest occupied band in the $\Gamma$ symmetry point to be aligned with 0 on the energy scale.

The above figure was plotted using Gnuplot and the following script:

#!/opt/local/bin/gnuplot

set terminal postscript enhanced eps color 16
set output 'band-gap.eps'

file1='rutile_example.restricted_band'
file2='rutile_example.dos'

set style line 1 lt 1 lw 3 pt 1 ps 1 lc rgb "#072c76"
set style line 2 lt 1 lw 3 pt 1 ps 1 lc rgb "#d61818"

set style line 100 lt 2 lw 1 pt 1 ps 1 lc rgb "#000000"

toev=27.211385
shift=0.233540E+00

set size 1.5,1.5
set multiplot

set size 1.0,1.5
set origin 0.0,0.0

set xrange [0:2.521406]

unset key

set ylabel "Energy [eV]"

set xtics ("{/Symbol G}" 0, "X" 0.361900, "M" 0.723800, \
"{/Symbol G}" 1.235603, "Z" 1.797606, "R" 2.159506, "A" 2.521406)

ymin=-20
ymax=10
set yrange [ymin:ymax]

set arrow from 0.361900,ymin to 0.361900,ymax nohead ls 100 # X
set arrow from 0.723800,ymin to 0.723800,ymax nohead ls 100 # M
set arrow from 1.235603,ymin to 1.235603,ymax nohead ls 100 # Gamma
set arrow from 1.797606,ymin to 1.797606,ymax nohead ls 100 # Z
set arrow from 2.159506,ymin to 2.159506,ymax nohead ls 100 # R

plot for [col=2:42] file1 using 1:((column(col)-shift)*toev) with lines ls 1

set size 0.5,1.5
set origin 1.0,0.0

set xtics 500
set xrange [0:1500]
unset ylabel
plot file2 u 2:((\$1-shift)*toev) with lines ls 2