import os

os.chdir("siesta_2")
import numpy as np
import sisl as si
import matplotlib.pyplot as plt

%matplotlib inline

Siesta — graphene

This tutorial will describe a complete walk-through of a large fraction of the sisl functionalities that may be related to the Siesta code.

Contrary to the ${\mathrm{H}}_2\mathrm{O}$ system this tutorial will emphasize the usefulness of performing bandstructures etc. directly in Python using sisl.

Creating the geometry

Our system of interest will be the smallest graphene cell. Instead of defining the atomic positions, the carbon atoms and supercell for graphene, we use a default implementation of graphene in sisl. There are a small selection of the typical geometries, including graphene.

graphene = si.geom.graphene(1.44)

Let us plot the geometry.

graphene.plot(axes="xy")

Now we need to create the input fdf file for Siesta:

open("RUN.fdf", "w").write(
    """%include STRUCT.fdf
SystemLabel siesta_2
PAO.BasisSize SZP
MeshCutoff 250. Ry
CDF.Save true
CDF.Compress 9
SaveHS true
SaveRho true
%block kgrid.MonkhorstPack
  61  1 1 0.
   1 61 1 0.
   0  0 1 0.
%endblock
"""
)
graphene.write("STRUCT.fdf")

Creating the electronic structure

Before proceeding, run Siesta to calculate the ground state electronic structure.

After having completed the Siesta run we may read the Hamiltonian to manipulate and extract different information. After reading the Hamiltonian it is obvious that a great deal of new data has been associated with the Hamiltonian.

fdf = si.get_sile("RUN.fdf")
H = fdf.read_hamiltonian()
print(H)

Hamiltonian{non-zero: 5994, orthogonal: False,
 Spin{unpolarized},
 Geometry{na: 2, no: 18,
  Atoms{species: 1,
   Atom{C, Z: 6, mass(au): 12.01000, maxR: 2.54650,
    AtomicOrbital{2sZ1, q0: 2.0, SphericalOrbital{l: 0, R: 2.175000000000011, q0: 2.0}},
    AtomicOrbital{2pyZ1, q0: 0.6666666666666666, SphericalOrbital{l: 1, R: 2.526900000000015, q0: 2.0}},
    AtomicOrbital{2pzZ1, q0: 0.6666666666666666, SphericalOrbital{l: 1, R: 2.526900000000015, q0: 2.0}},
    AtomicOrbital{2pxZ1, q0: 0.6666666666666666, SphericalOrbital{l: 1, R: 2.526900000000015, q0: 2.0}},
    AtomicOrbital{2dxyZ1P, q0: 0.0, SphericalOrbital{l: 2, R: 2.546500000000014, q0: 0.0}},
    AtomicOrbital{2dyzZ1P, q0: 0.0, SphericalOrbital{l: 2, R: 2.546500000000014, q0: 0.0}},
    AtomicOrbital{2dz2Z1P, q0: 0.0, SphericalOrbital{l: 2, R: 2.546500000000014, q0: 0.0}},
    AtomicOrbital{2dxzZ1P, q0: 0.0, SphericalOrbital{l: 2, R: 2.546500000000014, q0: 0.0}},
    AtomicOrbital{2dx2-y2Z1P, q0: 0.0, SphericalOrbital{l: 2, R: 2.546500000000014, q0: 0.0}}
   }: 2,
  },
  maxR: 2.54650,
  Lattice{nsc: [5 5 1],
   origin=[0.0000, 0.0000, 0.0000],
   A=[2.1600, 1.2471, 0.0000],
   B=[2.1600, -1.2471, 0.0000],
   C=[0.0000, 0.0000, 14.4000],
   bc=[Periodic,
       Periodic,
       Periodic]
  }
 }
}

Calculating DOS and PDOS

When we are dealing with periodic structures (such as graphene) it is imperative to calculate the density of states in a simple and efficient manner. Below we will calculate the DOS for a variety of Monkhorst-Pack grids to check the convergence of the DOS (it shouldn’t take more than a minute):

Warning

<p>
    This tutorial cannot fully converge the $k$-points. Please do the final convergence your self.
    You'll see little change after $\approx 60$.

E = np.linspace(-6, 4, 500)
for nk in [5, 11, 15, 21]:
    bz = si.MonkhorstPack(H, [nk, nk, 1])
    plt.plot(
        E,
        bz.apply.average.eigenvalue(wrap=lambda ev: ev.DOS(E)),
        label="nk={}".format(nk),
    )
plt.xlim(E[0], E[-1])
plt.ylim(0, None)
plt.xlabel(r"$E - E_F$ [eV]")
plt.ylabel(r"DOS [1/eV]")
plt.legend();

The default smearing method is the Gaussian smearing technique with $\sigma =0.1\mathrm{eV}$. Note that intrinsically the MonkhorstPack grid assumes time-reversal symmetry. I.e. $\mathbf{k}\equiv -\mathbf{k}$.

Now we may use the Monkhorst-Pack grid for 21 points (preferentially many more!) to find the projected DOS for some of the orbitals.

pdos_plot = H.plot.pdos(kgrid=[21, 21, 1], data_Erange=(-6, 4), Erange=[-6, 4], nE=500)

orb_groups = [
    {"l": 0, "name": "s", "color": "red"},
    {"l": 1, "m": [-1, 1], "name": "px + py", "color": "blue"},
    {"l": 1, "m": 0, "name": "pz", "color": "green"},
]
pdos_plot.update_inputs(groups=orb_groups)

As seen the $p_z$ orbitals are responsible for the DOS in a broad range of energies around the Fermi-level. This is one reason for the tight-binding models success with respect to graphene.

Another way to gain information is via the so-called fat-bands which basically is the PDOS scaling (no broadening) on each band for the quantities we are interested in. To plot the fat-bands we need the band-structure and a projection of each state onto the requested orbitals.

# Define the band-structure
bz = si.BandStructure(
    H,
    [[0] * 3, [2.0 / 3, 1.0 / 3, 0], [0.5, 0.5, 0], [1] * 3],
    400,
    names=[r"$\Gamma$", r"$K$", r"$M$", r"$\Gamma$"],
)

bz.plot.fatbands(groups=orb_groups, Erange=[-21, 10])

• all, black

• $p_z$, green

• $s$, red

• $p_x+p_y$, blue

Hamiltonian eigenstates

At this point we have plotted the $k$-averaged DOS, PDOS. We have also plotted the fat-bands (and thus the band-structure).

In addition to these things we can plot the real-space eigenstates. We first plot the $\mathrm{\Gamma }$-point for the first 2x2 unit-cell. This $k$-point has complete unit-cell periodicity and thus the plotted wavefunction should be fully periodic along all directions.

es = H.eigenstate()
idx_valence = (es.eig > 0).nonzero()[0][0] - 1
# Only select the valence band state
es = es.sub(idx_valence)

# Generate a grid encompassing a 2x2 graphene unit-cell
g = si.Grid(0.2, lattice=H.geometry.lattice.tile(2, 0).tile(2, 1))
# Calculate the real-space wavefunctions
es.wavefunction(g)

# Extract the wavefunction a few Ang above the graphene plane
# To do this we need to find the index of the corresponding z-plane.
# The Grid.index method is useful in this regard.
xyz = H.geometry.xyz[0, :].copy()
xyz[2] += 1.0
z_idx = g.index(xyz, axis=2)
x, y = np.mgrid[: g.shape[0], : g.shape[1]]
x, y = x * g.dcell[0, 0] + y * g.dcell[1, 0], x * g.dcell[0, 1] + y * g.dcell[1, 1]
plt.contourf(x, y, g.grid[:, :, z_idx])
xyz = H.geometry.tile(2, 0).tile(2, 1).xyz
plt.scatter(xyz[:, 0], xyz[:, 1], 20, c="k")
plt.xlabel(r"$x$ [Ang]")
plt.ylabel(r"$y$ [Ang]");

We will now try and plot the real-space wavefunction for a finite $k$. By choosing the $[1/2,0,0]$ point we know it must have a periodicity of 2 along the first lattice vector (this lattice vector is pointing right-up), and full periodicity along the second lattice vector. Since we have a finite $k$ the grid data-type must be complex because the eigenstates have complex components. And thus we will plot both the real and imaginary part.

_, axs = plt.subplots(1, 2, figsize=(16, 5))
es = H.eigenstate([1.0 / 2, 0, 0])
idx_valence = (es.eig > 0).nonzero()[0][0] - 1
es = es.sub(idx_valence)
g = si.Grid(0.2, dtype=np.complex128, lattice=H.geometry.lattice.tile(4, 0).tile(4, 1))
es.wavefunction(g)
x, y = np.mgrid[: g.shape[0], : g.shape[1]]
x, y = x * g.dcell[0, 0] + y * g.dcell[1, 0], x * g.dcell[0, 1] + y * g.dcell[1, 1]
axs[0].contourf(x, y, g.grid[:, :, z_idx].real)
axs[1].contourf(x, y, g.grid[:, :, z_idx].imag)
xyz = H.geometry.tile(4, 0).tile(4, 1).xyz
for ax in axs:
    ax.scatter(xyz[:, 0], xyz[:, 1], 20, c="k")
    ax.set_xlabel(r"$x$ [Ang]")
    ax.set_ylabel(r"$y$ [Ang]");

info:0: SislInfo:

wavefunction: k != Gamma is currently untested!