### Introduction

For this exercise we will try to reproduce a band structure and density of states of TiO_{2}
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 (TiO_{2}) 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
task band jobType
```

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 TiO_{2} in the rutile structure. We will specify
the dimensions of the unit cell in Angstroms and angles in degrees. TiO_{2} 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
TiO_{2} 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
task band energy
```

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 TiO_{2} 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
task band energy
nwpw
virtual 40
brillouin_zone
zone_name tetragonalzone
path simple-tetragonal gamma x m gamma z r a
end
zone_structure_name tetragonalzone
end
task band structure
```

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
task band dos
```

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
task band energy
set dos:npoints 1000
nwpw
virtual 40
dos-grid 11 11 11
end
task band dos
```

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 TiO_{2} 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
```