In this notebook we will see, how to calculate band structures using pyiron.
from pyiron import Project
from ase.spacegroup import crystal
import matplotlib.pyplot as plt
import seekpath as sp
import numpy as np
pr = Project("band_structure")
First we see how seekpath works! Therefore we create firts a structure using pyiron.
structure_Fe = pr.create_structure("Fe", "bcc", 2.81)
structure_Fe.plot3d()
_ColormakerRegistry()
NGLWidget()
structure_Fe
Fe: [0. 0. 0.] Fe: [1.405 1.405 1.405] pbc: [ True True True] cell: [[2.81000000e+00 0.00000000e+00 0.00000000e+00] [1.72062875e-16 2.81000000e+00 0.00000000e+00] [1.72062875e-16 1.72062875e-16 2.81000000e+00]]
For seekpath we need a tuple containing
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': [[1.6776982741209693e-16, 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.
We use the same structure as before!
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()
NGLWidget()
structure_Fe_sp
Fe: [0. 0. 0.] pbc: [ True True True] 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')]}
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.
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.load_default_groups()
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)
ham_spx_chg.run()
The job Fe_spx_CHG was saved and received the ID: 4
pr.get_jobs_status()
finished 1 Name: status, dtype: int64
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")
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/parameters.py:285: UserWarning: The input in GenericParameters changed, while the state of the job was already finished. "The input in GenericParameters changed, while the state of the job was already finished."
A parameter usefull for BS calculations. Look at the sphinx manual for details.
ham_spx_bs.input["nSloppy"] = 6
ham_spx_bs.run()
The job Fe_spx_BS was saved and received the ID: 5
pr.get_jobs_status()
finished 2 Name: status, dtype: int64
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
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);