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
%config InlineBackend.figure_format = 'retina'
pr = Project("band_structure")
pr.remove_jobs(recursive=True)
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': [[4.160166955977647e-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: 2
--------------------------------------------------------------------------- AssertionError Traceback (most recent call last) <ipython-input-18-694314fe4b23> in <module> ----> 1 ham_spx_chg.run() ~/Software/pyiron/pyiron/base/job/generic.py in run(self, run_again, repair, debug, run_mode) 662 status = "created" 663 if status == "initialized": --> 664 self._run_if_new(debug=debug) 665 elif status == "created": 666 que_id = self._run_if_created() ~/Software/pyiron/pyiron/base/job/generic.py in _run_if_new(self, debug) 1305 else: 1306 self._create_job_structure(debug=debug) -> 1307 self.run() 1308 1309 def _run_if_created(self): ~/Software/pyiron/pyiron/base/job/generic.py in run(self, run_again, repair, debug, run_mode) 664 self._run_if_new(debug=debug) 665 elif status == "created": --> 666 que_id = self._run_if_created() 667 if que_id: 668 self._logger.info( ~/Software/pyiron/pyiron/base/job/generic.py in _run_if_created(self) 1322 self.run_if_manually() 1323 elif self.server.run_mode.modal: -> 1324 self.run_static() 1325 elif self.server.run_mode.non_modal or self.server.run_mode.thread: 1326 self.run_if_non_modal() ~/Software/pyiron/pyiron/base/job/generic.py in run_static(self) 764 "{}, status: {}, output: {}".format(self.job_info_str, self.status, out) 765 ) --> 766 self.run() 767 if job_crashed: 768 self.status.aborted = True ~/Software/pyiron/pyiron/base/job/generic.py in run(self, run_again, repair, debug, run_mode) 677 self._run_if_running() 678 elif status == "collect": --> 679 self._run_if_collect() 680 elif status == "suspend": 681 self._run_if_suspended() ~/Software/pyiron/pyiron/base/job/generic.py in _run_if_collect(self) 1393 status is set to 'finished' 1394 """ -> 1395 self.collect_output() 1396 self.collect_logfiles() 1397 if self.job_id is not None: ~/Software/pyiron/pyiron/sphinx/base.py in collect_output(self, force_update) 1274 Collects the outputs and stores them to the hdf file 1275 """ -> 1276 self._output_parser.collect(directory=self.working_directory) 1277 self._output_parser.to_hdf(self._hdf5, force_update=force_update) 1278 ~/Software/pyiron/pyiron/sphinx/base.py in collect(self, directory) 2305 self.collect_energy_struct(file_name="energy-structOpt.dat", 2306 cwd=directory) -> 2307 self.collect_sphinx_log(file_name="sphinx.log", cwd=directory) 2308 self.collect_relaxed_hist(file_name="relaxHist.sx", cwd=directory) 2309 self._job.compress() ~/Software/pyiron/pyiron/sphinx/base.py in collect_sphinx_log(self, file_name, cwd, check_consistency) 2054 self._job.status.aborted = True 2055 raise AssertionError("SPHInX did not enter the main loop; \ -> 2056 output not collected") 2057 if not np.any(["Program exited normally." in line 2058 for line in log_file]): AssertionError: SPHInX did not enter the main loop; output not collected
pr.get_jobs_status()
aborted 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)
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: 3
--------------------------------------------------------------------------- AssertionError Traceback (most recent call last) <ipython-input-23-31eb5b79eaf6> in <module> ----> 1 ham_spx_bs.run() ~/Software/pyiron/pyiron/base/job/generic.py in run(self, run_again, repair, debug, run_mode) 662 status = "created" 663 if status == "initialized": --> 664 self._run_if_new(debug=debug) 665 elif status == "created": 666 que_id = self._run_if_created() ~/Software/pyiron/pyiron/base/job/generic.py in _run_if_new(self, debug) 1305 else: 1306 self._create_job_structure(debug=debug) -> 1307 self.run() 1308 1309 def _run_if_created(self): ~/Software/pyiron/pyiron/base/job/generic.py in run(self, run_again, repair, debug, run_mode) 664 self._run_if_new(debug=debug) 665 elif status == "created": --> 666 que_id = self._run_if_created() 667 if que_id: 668 self._logger.info( ~/Software/pyiron/pyiron/base/job/generic.py in _run_if_created(self) 1322 self.run_if_manually() 1323 elif self.server.run_mode.modal: -> 1324 self.run_static() 1325 elif self.server.run_mode.non_modal or self.server.run_mode.thread: 1326 self.run_if_non_modal() ~/Software/pyiron/pyiron/base/job/generic.py in run_static(self) 764 "{}, status: {}, output: {}".format(self.job_info_str, self.status, out) 765 ) --> 766 self.run() 767 if job_crashed: 768 self.status.aborted = True ~/Software/pyiron/pyiron/base/job/generic.py in run(self, run_again, repair, debug, run_mode) 677 self._run_if_running() 678 elif status == "collect": --> 679 self._run_if_collect() 680 elif status == "suspend": 681 self._run_if_suspended() ~/Software/pyiron/pyiron/base/job/generic.py in _run_if_collect(self) 1393 status is set to 'finished' 1394 """ -> 1395 self.collect_output() 1396 self.collect_logfiles() 1397 if self.job_id is not None: ~/Software/pyiron/pyiron/sphinx/base.py in collect_output(self, force_update) 1274 Collects the outputs and stores them to the hdf file 1275 """ -> 1276 self._output_parser.collect(directory=self.working_directory) 1277 self._output_parser.to_hdf(self._hdf5, force_update=force_update) 1278 ~/Software/pyiron/pyiron/sphinx/base.py in collect(self, directory) 2305 self.collect_energy_struct(file_name="energy-structOpt.dat", 2306 cwd=directory) -> 2307 self.collect_sphinx_log(file_name="sphinx.log", cwd=directory) 2308 self.collect_relaxed_hist(file_name="relaxHist.sx", cwd=directory) 2309 self._job.compress() ~/Software/pyiron/pyiron/sphinx/base.py in collect_sphinx_log(self, file_name, cwd, check_consistency) 2054 self._job.status.aborted = True 2055 raise AssertionError("SPHInX did not enter the main loop; \ -> 2056 output not collected") 2057 if not np.any(["Program exited normally." in line 2058 for line in log_file]): AssertionError: SPHInX did not enter the main loop; output not collected
pr.get_jobs_status()
aborted 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);