import os

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

%matplotlib inline

info:0: SislInfo: Please install tqdm (pip install tqdm) for better looking progress bars

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 = geom.graphene(1.44)

Let us plot the geometry.

graphene.plot(axes="xy")

_ipython_canary_method_should_not_exist_ False
_ipython_display_ True

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 = get_sile("RUN.fdf")
H = fdf.read_hamiltonian()
print(H)

Hamiltonian{non-zero: 5994, orthogonal: False,
 Spin{unpolarized, kind=f},
 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):

E = np.linspace(-6, 4, 500)
for nk in [21, 41, 61, 81]:
    bz = 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();

We see that the k-points are converged for approximately 61 k-points using the default smearing method. 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 81 points to find the projected DOS for some of the orbitals.

pdos_plot = H.plot.pdos(kgrid=[81, 81, 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)

_ipython_canary_method_should_not_exist_ False
_ipython_display_ True

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 = 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])

_ipython_canary_method_should_not_exist_ False
_ipython_display_ True

• 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 \(\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 = 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 = 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!