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

  1. The cell in \(3\times3\) array

  2. The scaled positions

  3. 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=100and 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);
../_images/sol_4_bandstructure_55_0.png