# Exercise 4: Band structure calculations¶

In this notebook we will see, how to calculate band structures using pyiron.

```
from pyiron.project import Project
from ase.spacegroup import crystal
import matplotlib.pyplot as plt
import seekpath as sp
import numpy as np
```

```
pr = Project("demo/band_structure")
```

```
# pr.remove_jobs(recursive=True)
```

## Structures with seekpath¶

First we see how seekpath works! Therefore we create firts a structure using pyiron.

## Create structure with pyiron¶

```
structure_Fe = pr.create_structure("Fe", "bcc", 2.81)
```

```
structure_Fe.plot3d()
```

```
structure_Fe
```

```
Fe: [0. 0. 0.]
Fe: [1.405 1.405 1.405]
pbc: [ True True True]
cell:
Cell([[2.81, 0.0, 0.0], [1.7206287528020313e-16, 2.81, 0.0], [1.7206287528020313e-16, 1.7206287528020313e-16, 2.81]])
```

## Create structure with seekpath¶

For seekpath we need a tuple containing

The cell in \(3\times3\) array

The scaled positions

List of

`ints`

to distinguish the atom types (indices of pyiron structure) as input structure.

```
input_sp = (structure_Fe.cell, structure_Fe.get_scaled_positions(), structure_Fe.indices)
```

Just to see how the output looks like, let us do…

```
sp.get_path(input_sp)
```

```
{'point_coords': {'GAMMA': [0.0, 0.0, 0.0],
'H': [0.5, -0.5, 0.5],
'P': [0.25, 0.25, 0.25],
'N': [0.0, 0.0, 0.5]},
'path': [('GAMMA', 'H'),
('H', 'N'),
('N', 'GAMMA'),
('GAMMA', 'P'),
('P', 'H'),
('P', 'N')],
'has_inversion_symmetry': True,
'augmented_path': False,
'bravais_lattice': 'cI',
'bravais_lattice_extended': 'cI1',
'conv_lattice': array([[2.81, 0. , 0. ],
[0. , 2.81, 0. ],
[0. , 0. , 2.81]]),
'conv_positions': array([[0. , 0. , 0. ],
[0.5, 0.5, 0.5]]),
'conv_types': array([0, 0], dtype=int32),
'primitive_lattice': array([[-1.405, 1.405, 1.405],
[ 1.405, -1.405, 1.405],
[ 1.405, 1.405, -1.405]]),
'primitive_positions': array([[0., 0., 0.]]),
'primitive_types': array([0], dtype=int32),
'reciprocal_primitive_lattice': [[-0.0, 2.236009006113732, 2.236009006113732],
[2.236009006113732, 0.0, 2.236009006113732],
[2.236009006113732, 2.236009006113732, 0.0]],
'inverse_primitive_transformation_matrix': array([[0, 1, 1],
[1, 0, 1],
[1, 1, 0]]),
'primitive_transformation_matrix': array([[-0.5, 0.5, 0.5],
[ 0.5, -0.5, 0.5],
[ 0.5, 0.5, -0.5]]),
'volume_original_wrt_conv': 0.9999999999999999,
'volume_original_wrt_prim': 1.9999999999999998,
'spacegroup_number': 229,
'spacegroup_international': 'Im-3m'}
```

The code creates automatically the conventional and primitive cell with all high-symmetry points and a suggested path taking all high-symmetry paths into account.

**Keep in mind:** The high-symmetry points and paths make only sence for the primitive cell! Therefore we run all calculations in the primitive cell created by seekpath.

## Create a structure¶

We use the same structure as before!

## Create new structure (primitive cell) with high-symmetry points and paths¶

For the following command all arguments valid for seekpath are supported. Look at the docstring and at seekpath.

```
structure_Fe_sp = structure_Fe.create_line_mode_structure()
```

```
structure_Fe_sp.plot3d()
```

```
structure_Fe_sp
```

```
Fe: [0. 0. 0.]
pbc: [ True True True]
cell:
Cell([[-1.405, 1.405, 1.405], [1.405, -1.405, 1.405], [1.405, 1.405, -1.405]])
```

We see, that the structure is now the primitive cell containing only one atom.

```
structure_Fe_sp.get_high_symmetry_points()
```

```
{'GAMMA': [0.0, 0.0, 0.0],
'H': [0.5, -0.5, 0.5],
'P': [0.25, 0.25, 0.25],
'N': [0.0, 0.0, 0.5]}
```

```
structure_Fe_sp.get_high_symmetry_path()
```

```
{'full': [('GAMMA', 'H'),
('H', 'N'),
('N', 'GAMMA'),
('GAMMA', 'P'),
('P', 'H'),
('P', 'N')]}
```

The path is stored like this. Here you can also add paths to the dictionary.

Each tuple gives a start and end point for this specific trace. Thus also disonnected paths are possible to calculate.

```
structure_Fe_sp.add_high_symmetry_path({"my_path": [("GAMMA", "H"), ("GAMMA", "P")]})
```

```
structure_Fe_sp.get_high_symmetry_path()
```

```
{'full': [('GAMMA', 'H'),
('H', 'N'),
('N', 'GAMMA'),
('GAMMA', 'P'),
('P', 'H'),
('P', 'N')],
'my_path': [('GAMMA', 'H'), ('GAMMA', 'P')]}
```

## Create jobs¶

We need two jobs for a band structure! The first gives us the correct Fermi energy and the charge densities used for the second calculations.

## Create job for charge density¶

This is only a small example for BS calculations. Could be that the input parameter like cutoff etc. does not make much sense… for real physics…

```
def setup_hamiltonian_sphinx(project, jobname, structure, chgcar_file=""):
#version 1.0 (08.03.2019)
#Name und typ
ham = project.create_job(job_type='Sphinx', job_name=jobname)
#parameter für xc functional
ham.exchange_correlation_functional = 'PBE'
#struktur
ham.structure = structure
ham.set_encut(450)
ham.set_empty_states(6)
ham.set_convergence_precision(electronic_energy=1e-8)
ham.set_occupancy_smearing(width=0.2)
#parameter für kpoints
ham.set_kpoints([8, 8, 8])
return ham
```

```
ham_spx_chg = setup_hamiltonian_sphinx(pr, "Fe_spx_CHG", structure_Fe_sp)
```

### Run it!¶

```
ham_spx_chg.run()
```

```
The job Fe_spx_CHG was saved and received the ID: 127
```

```
pr.get_jobs_status()
```

```
finished 1
Name: status, dtype: int64
```

## Create second job¶

We restart the fist job with the following command. Then the charge density of the first job is taken for the second!

```
ham_spx_bs = ham_spx_chg.restart_for_band_structure_calculations(job_name="Fe_spx_BS")
```

### Set line mode for k-points¶

To set the correct path, we have to give the name of the path (in our example either `full`

or `my_path`

) and the number of points for each subpath (would be for `n_path=100`

and `path_name="my_path"`

200 k-points in total)

```
ham_spx_bs.set_kpoints(scheme="Line", path_name="full", n_path=100)
```

```
/srv/conda/envs/notebook/lib/python3.7/site-packages/pyiron_base/generic/inputlist.py:325: UserWarning: The input in Group changed, while the state of the job was already finished.
"finished.".format(cls.__name__)
```

A parameter usefull for BS calculations. Look at the sphinx manual for details.

```
ham_spx_bs.input["nSloppy"] = 6
```

### Run it!¶

```
ham_spx_bs.run()
```

```
The job Fe_spx_BS was saved and received the ID: 135
```

```
pr.get_jobs_status()
```

```
finished 2
Name: status, dtype: int64
```

## Store the data!¶

The energy values are stored in the following paths of the hdf5 file.

```
energy_sphinx = ham_spx_bs['output/generic/dft/bands_eigen_values'][-1]
ef_sphinx = ham_spx_chg['output/generic/dft/bands_e_fermi'][-1]
energy_sphinx -= ef_sphinx
```

## Plot it!¶

Now we can easily plot it!

```
plt.plot(energy_sphinx[:-1], 'b-')
plt.axhline(y=0, ls='--', c='k')
plt.xlim(0, len(energy_sphinx));
plt.ylim(-10,40);
```