[1]:
import numpy as np
import sisl as si
from sisl.viz import merge_plots
from sisl.viz.processors.math import normalize
import matplotlib.pyplot as plt
%matplotlib inline
Electronic structure calculation – a walk-through
This tutorial will describe a complete walk-through of a large fraction of the sisl
functionalities. It will show you how to generated default geometries, constructing Hamiltonians, calculating eigenstates and plotting various physical quantities.
Creating the geometry to investigate
Our system of interest will be graphene. Instead of creating a graphene flake, or the primary unit-cell of graphene, we will create a vacancy in graphene. We will start by creating a graphene flake
[2]:
graphene = si.geom.graphene().tile(6, 0).tile(6, 1)
This does a lot of things behind the scenes:
si.geom.graphene
:create atomic coordinates of pristine graphene with a default bond-length of \(1.42\,\mathrm{Å}\).
create pristine graphene unit cell, by default this will create a supercell with a size
3x3
, i.e. a nearest neighbour unit-cell.assign a carbon atom with a default of one orbital per atom as the basis
Geometry.tile
tiles the geometry(reps, axis)
byreps
times along the unit cell axisaxis
By printing the object one gets basic information regarding the geometry, such as 1) number and species of atoms, 2) number of orbitals, 3) orbitals associated with each atom and 4) number of supercells.
[3]:
print(graphene)
Geometry{na: 72, no: 72,
Atoms{species: 1,
Atom{C, Z: 6, mass(au): 12.01070, maxR: 1.43420,
Orbital{R: 1.43420, q0: 0.0}
}: 72,
},
maxR: 1.43420,
Lattice{nsc: [3 3 1],
origin=[0.0000, 0.0000, 0.0000],
A=[12.7800, -7.3785, 0.0000],
B=[12.7800, 7.3785, 0.0000],
C=[0.0000, 0.0000, 20.0000],
bc=[Periodic,
Periodic,
Unknown]
}
}
In this example we have na=72
atoms, each have one orbital, hence the total number of orbitals is no=72
. The description of the atomic specie in the geometry tells us we have a carbon atom, with a single orbital with a radius of \(1.4342\,\mathrm{Å}\). The number of supercells are [3, 3, 1]
which means cells {-1, 0, 1}
along the first and second lattice are taken into account.
Later we will look into the details of orbitals associated with atoms and how they may be used for wavefunctions etc.
Lets visualize the atomic positions (here adding atomic indices)
[4]:
graphene.plot(axes="xy")
Removing an atom can be done with Geometry.remove
. The routine takes an index, or a list of indices of the atoms to be removed. For instance removing the first atom will result in the following geometry (red atom is the removed atom)
[5]:
coord = graphene.sub(0)
coord.plot(axes="xy", atoms_style={"color": "red"}).merge(
graphene.remove(0).plot(axes="xy")
)
For the following it doesn’t matter which atom we remove (since it is peridiodic), however for visualization purposes we will remove an atom in the middle of the unit cell.
Using sisl
it is easy to find atoms close to specific positions. The middle of the atomic coordinates is also the center of atomic coordinates, here Geometry.center
is useful. The complementary method Geometry.close
finds all atomic indices close to a given position or atom. By default, Geometry.close
determines all atoms within a radius equal to the maximum orbital radius. Here we explicitly set the search radius to \(1.5\,\mathrm{Å}\).
[6]:
xyz_center = graphene.center(what="xyz")
indices = graphene.close(xyz_center, 1.5)
index = indices[0]
system = graphene.remove(index)
graphene.plot(
axes="xy",
atoms_style=[
{"opacity": 0.5}, # Default style for all atoms
{
"atoms": indices,
"color": "black",
"size": 1.2,
"opacity": 1,
}, # Styling for indices_close_to_center on top of defaults.
{
"atoms": index,
"color": "red",
"size": 1,
"opacity": 1,
}, # Styling for center_atom_index on top of defaults.
],
)
Creating the electronic structure
To calculate eigenstates, DOS and other physical quantities from the Hamiltonian we need to setup the Hamiltonian. Tthis is done by passing a Geometry
to a Hamiltonian
class
[7]:
H = si.Hamiltonian(system)
print(H)
Hamiltonian{non-zero: 0, orthogonal: True,
Spin{unpolarized, kind=f},
Geometry{na: 71, no: 71,
Atoms{species: 1,
Atom{C, Z: 6, mass(au): 12.01070, maxR: 1.43420,
Orbital{R: 1.43420, q0: 0.0}
}: 71,
},
maxR: 1.43420,
Lattice{nsc: [3 3 1],
origin=[0.0000, 0.0000, 0.0000],
A=[12.7800, -7.3785, 0.0000],
B=[12.7800, 7.3785, 0.0000],
C=[0.0000, 0.0000, 20.0000],
bc=[Periodic,
Periodic,
Unknown]
}
}
}
Geometry
information it informs us that it is an orthogonal basis (sisl
also allows non-orthogonal basis’). The spin-configuration is an unpolarized configuration (see Spin
for details).non-zero = 0
specifies that there are no associated Hamiltonian elements stored in this Hamiltonian object.Hamiltonian.construct
lets one specify the Hamiltonian elements in a consistent way. However, note that the Hamiltonian
objet may be used as though it was a matrix, i.e. Hamiltonian[0, 1] = a
will set the hopping element from the 0th orbital to the 1st orbital to a
.construct
method is a list of radii and an accompanying list of energies. Here r
tells sisl
to find all atoms within a sphere of \(0.1\,\mathrm{Å}\) from each atom and set the corresponding element to \(0\), secondly all atoms within a spherical radius of \(0.1\,\mathrm{Å}\) to \(1.44\,\mathrm{Å}\) are given the matrix element \(-2.7\,\mathrm{eV}\)[8]:
r = (0.1, 1.44)
t = (0.0, -2.7)
H.construct([r, t])
print(H)
Hamiltonian{non-zero: 281, orthogonal: True,
Spin{unpolarized, kind=f},
Geometry{na: 71, no: 71,
Atoms{species: 1,
Atom{C, Z: 6, mass(au): 12.01070, maxR: 1.43420,
Orbital{R: 1.43420, q0: 0.0}
}: 71,
},
maxR: 1.43420,
Lattice{nsc: [3 3 1],
origin=[0.0000, 0.0000, 0.0000],
A=[12.7800, -7.3785, 0.0000],
B=[12.7800, 7.3785, 0.0000],
C=[0.0000, 0.0000, 20.0000],
bc=[Periodic,
Periodic,
Unknown]
}
}
}
Now the Hamiltonian has \(281=71\cdot4-3\) non-zero elements (as expected).
Hamiltonian eigenstates
To begin with we calculate the \(\Gamma\)-point eigenstates and plot a subset of the eigenstates’ norm on the geometry.
[9]:
es = H.eigenstate()
# Reduce the contained eigenstates to only 3 states around the Fermi-level
es_fermi = es.sub(range(len(H) // 2 - 1, len(H) // 2 + 2))
plots = [
system.plot(axes="xy", atoms_style=[{"size": n * 20, "color": c}])
for n, c in zip(es_fermi.norm2(projection="orbital"), ("red", "blue", "green"))
]
merge_plots(*plots, composite_method="subplots", cols=3)
The Hamiltonian.eigenstate
(with an optional \(k\)-point argument) routine returns an EigenstateElectron
object which holds the eigenvalues and eigenvectors for a given \(k\)-point. This object can perform several advanced calculations:
EigenstateElectron.DOS
: calculate the DOS at a given set of energy values (in eV), additionally one can pass a distribution function if the default Gaussian with \(\sigma=0.1\,\mathrm{eV}\) is not wanted.EigenstateElectron.PDOS
: calculate the projected DOS at a given set of energy values (in eV), additionally one can pass a distribution function if the default Gaussian with \(\sigma=0.1\,\mathrm{eV}\) is not wanted.EigenstateElectron.wavefunction
: add all contained eigenstates to a passed real-space grid.
Lets first try and calculate the DOS for a given set of energies in the range \(-4\,\mathrm{eV}\) : \(\mathrm{4}\,\mathrm{eV}\)
[10]:
E = np.linspace(-4, 4, 400)
plt.plot(E, es.DOS(E))
plt.xlabel(r"$E - E_F$ [eV]")
plt.ylabel(r"DOS at $\Gamma$ [1/eV]");
The projected DOS (in this case) can aptly be plotted on individual atoms as will be seen in the following. We will integrate the PDOS in the range \(-1\,\mathrm{eV}\) to \(-0.5\,\mathrm{eV}\)
[11]:
E = np.linspace(-1, -0.5, 100)
dE = E[1] - E[0]
PDOS = es.PDOS(E).sum((0, 2)) * dE # perform integration
system.plot(axes="xy", atoms_style={"size": normalize(PDOS, 0, 1)})
# plt.scatter(system.xyz[:, 0], system.xyz[:, 1], 500 * PDOS);
# plt.scatter(xyz_remove[0], xyz_remove[1], c='k', marker='*'); # mark the removed atom
This highlights a somewhat localized state around the missing atom.
Brillouin-zone calculations
The above DOS and PDOS analysis are useful for investigating a single \(k\)-point at a time. However, they are incomplete in the sense of the full Brillouin zone. To leverage this sisl
, implements several classes to handle Brillouin-zones.
In the following we will show how to perform band structure calculations as well as performing \(k\)-averaged quantities in the Brillouin zone.
Bandstructure
sisl
calculating the band-structure is as easy as anything else.[12]:
band = si.BandStructure(
H,
[[0, 0, 0], [0, 0.5, 0], [1 / 3, 2 / 3, 0], [0, 0, 0]],
400,
[r"Gamma", r"M", r"K", r"Gamma"],
)
Note that the integer 400
determines the total number of \(k\)-points on the full band. One can define an explicit number between the different points, however we highly encourage only specifying one integer as the divisions will be determined based on the length in the reciprocal space between the points and thus the physical distance in the Brillouin zone will be correct.
A word of note on the BrillouinZone
objects.
The BrillouinZone
objects are extremely handy because they allow to directly call any routine inherent to the passed object. If calling routine band.a()
it is equivalent to:
for k in band: yield band.parent.a(k=k)
Note that BrillouinZone
defaults to use BrillouinZone.apply.iter
.
However, for large systems this may result in memory problems in which case it may be necessary to return an iterator. To circumvent this we can tell the Brillouin zone object to (always) return a list instead:
b_list = band.apply.list
as_list = b_list.a()
Another option is to return \(k\)-point averaged quantities which is roughly equivalent to (internally it is done with minimal memory usage):
band.apply.average.a() == sum([band.parent.a(k=k) * band.weight[i] for i, k in enumerate(band)])
Now we can calculate the band structure. A band-structure requires all eigenvalues and thus we ask the BrillouinZone
object to return all values using apply.array
.
[13]:
band.plot(Erange=[-3, 3])
Calculating \(k\)-averaged quantities
Now we are in a position to calculate a subset of quantities from the Hamiltonian. Before proceeding we will note that the Hamiltonian also implements the DOS
and PDOS
methods (equivalent to Hamiltonian.eigenvalue().DOS()
), hence to calculate these as \(k\)-averaged quantities we can create a Brillouin zone object with proper weights, say a Monkhorst-Pack grid, and calculate the averaged quantities.
[14]:
bz = si.MonkhorstPack(H, [35, 35, 1])
bz_average = (
bz.apply.average
); # specify the Brillouin zone to perform an average of subsequent method calls
[15]:
E = np.linspace(-4, 4, 1000)
plt.plot(E, bz_average.eigenstate(wrap=lambda es: es.DOS(E)))
plt.xlabel("$E - E_F$ [eV]")
plt.ylabel("DOS [1/eV]");
We can also plot the projected DOS integrated around the Fermi level on the atoms
[16]:
E = np.linspace(-1, 1, 1000)
dE = E[1] - E[0]
PDOS = bz_average.eigenstate(wrap=lambda es: es.PDOS(E).sum((0, 2))) * dE
[17]:
system.plot(axes="xy", atoms_style={"size": normalize(PDOS, 0, 1)})
Plotting eigenstates on a real space grid
sisl
also implements methods to plot orbitals on grids. Often it may also be advantageous to plot simple orbitals to check their appearence. sisl
implements a simple variation of spherical atomic orbitals. Other orbitals may easily be added, if so desired.
Since we require orbitals to be zero at some maximum cut-off \(r_\mathrm{max}\) a radial function is used to cutoff the spherical harmonics. In this case we will simply use an exponential function with a cutoff radius of \(1.6\,\mathrm{Å}\)
[18]:
r = np.linspace(0, 1.6, 700)
f = 5 * np.exp(-r * 5)
print("Normalization: {}".format(f.sum() * (r[1] - r[0])))
plt.plot(r, f)
plt.ylim([0, None])
orb = si.SphericalOrbital(1, (r, f));
Normalization: 1.0053998295349618
To calculate the orbital on a 3D grid, the Orbital.toGrid
function is available. Here we create an orbital with azimuthal quantum number \(l=1\) and by default it has angular moment \(0\), i.e. the \(2p_z\) orbital. To create a \(2p_x\) or \(2p_y\) orbital requires the use of an AtomicOrbital
to denote the \(m\) quantum number. Below we create the grid that describes the \(p_z\) orbital and plot the \(yz\) plane at \(x=0\) (with respect to the position
of the orbital).
Note that the real-space orbital will not be normalized and thus \(\int|\psi(\mathbf r)|^2\neq 1\), as it should. The following analysis is thus a qualitative analysis.
[19]:
grid = orb.toGrid()
index = grid.index(-grid.origin)
[20]:
plt.imshow(grid.grid[index[0], :, :].T);
This simple orbital model may be used to plot the real-space tight-binding eigenstates. Our first task will be to tell the geometry that the orbital associated with the Carbon atoms is the \(2p_z\) orbital as noted above. First we create the Carbon atom, then we replace the atom in the system geometry.
[21]:
C = si.Atom(6, orb)
print(system.atoms)
system.atoms.replace(system.atoms[0], C)
print(system.atoms)
Atoms{species: 1,
Atom{C, Z: 6, mass(au): 12.01070, maxR: 1.43420,
Orbital{R: 1.43420, q0: 0.0}
}: 71,
}
Atoms{species: 1,
Atom{C, Z: 6, mass(au): 12.01070, maxR: 1.55180,
SphericalOrbital{l: 1, R: 1.5517999999999998, q0: 0.0}
}: 71,
}
Note the difference between the two print(system)
statements, before the replacement, an intrinsic Orbital
object describes the orbital, in effect no knowledge other than the radius. After replacement, the spherical orbital with azimuthal angular moment \(l=1\) is replaced. Now we can plot the real-space grid for an eigenstate.
Lets plot one of the \(\Gamma\)-point eigenstates close to the Fermi-level.
[22]:
es = H.eigenstate(dtype=np.float64).sub([len(H) // 2 + 1])
grid = si.Grid(0.075, lattice=H.lattice)
es.wavefunction(grid)
[23]:
index = grid.index([0, 0, 0.1])
grid = grid.sub(index, 2)
[24]:
grid.plot(axes="xy")
[25]:
grid.plot(axes="ab")