Setting up the primitive cell

In this tutorial, we show how to set up the primitive cell from the PrimitiveCell class taking monolayer graphene as the example. We will discuss the geometric properties of graphene first, then set up the model step-by-step. The script in this tutorial can be found at examples/prim_cell/model/graphene_diamond.py.

Geometric properties

Monolayer graphene has a hexagonal lattice, with lattice paprameters \(a=b=2.46 \overset{\circ}{\mathrm {A}}\) and \(\alpha=\beta=90^\circ\). The angle between lattice vectors \(a_1\) and \(a_2\), namely \(\gamma\), can be either \(60^\circ\) or \(120^\circ\), depending on the choice of lattice vectors, as shown in the figure:

../../_images/graph_lattice.png

Crystalline structure of monolayer graphene. Grey spheres indicate the carbon atoms. \((a_1, a_2)\) and \((a_1, a_2\prime)\) are two equivalent sets of lattice vectors, with the angle \(\gamma\) being either \(60^\circ\) or \(120^\circ\). The corresponding primitive cells are indicated with red and blue diamonds, respectively.

Each primitive cell of monolayer graphene has two carbon atoms located at

  • \(\tau_1 = \frac{1}{3}a_1 + \frac{1}{3}a_2\)

  • \(\tau_2 = \frac{2}{3}a_1 + \frac{2}{3}a_2\)

if \(\gamma=60^\circ\) and

  • \(\tau_1 = \frac{2}{3}a_1 + \frac{1}{3}a_2\)

  • \(\tau_2 = \frac{1}{3}a_1 + \frac{2}{3}a_2\)

if \(\gamma=120^\circ\). We call the coordinates of \((\frac{1}{3}, \frac{1}{3}) (\frac{2}{3}, \frac{2}{3})\) and \((\frac{2}{3}, \frac{1}{3}) (\frac{1}{3}, \frac{2}{3})\) as fractional coordinates. If you are familiar with density functional theory programs, e.g. VASP and Quantum ESPRESSO, you will immediately see that they are the direct or crystal coordinates. The advantage of fractional coordinates is that they are compact and independent of the orientation of lattice vectors. The PrimitiveCell class of TBPLaS uses fractional coordinates internally, yet there are utilities to convert them from/to Cartesian coordinates.

Create empty cell

First of all, we need to import all necessary packages:

import math
import numpy as np
import tbplas as tb

To create an empty primitive cell, we need to define of lattice vectors. We choose \(\gamma=60^\circ\). Since all cells are implemented as 3-dimensional within TBPLaS, we need to specify an arbitrary lattice constant c, which is taken to be \(10\overset{\circ}{\mathrm {A}}\). The lattice vectors can be generated by:

vectors = tb.gen_lattice_vectors(a=2.46, b=2.46, c=10.0, gamma=60)

Alternatively, you can input the coordinates manually:

a = 2.46
c = 10.0
a_half = a * 0.5
sqrt3 = math.sqrt(3)

vectors = np.array([
    [a, 0, 0,],
    [a_half, sqrt3*a_half, 0],
    [0, 0, c]
])

Note that the orientation of the lattice vectors is quite arbitrary, and we have \(a_1\) parallel to \(x\)-axis just for convenience. We can also make \(a_1\) and \(a_2\) to be symmetric about \(x\)-axis.

vectors = np.array([
    [sqrt3*a_half, -a_half, 0],
    [sqrt3*a_half, a_half, 0],
    [0, 0, c]
])

or symmetric about \(y\)-axis:

vectors = np.array([
    [a_half, sqrt3*a_half, 0],
    [-a_half, sqrt3*a_half, 0],
    [0, 0, c]
])

When the lattice vectors are ready, we can create an empty primitive cell by:

cell = tb.PrimitiveCell(vectors)

Add orbitals

Since we choose \(\gamma=60^\circ\), the two carbon atoms are then located at \((\frac{1}{3}, \frac{1}{3})\) and \((\frac{2}{3}, \frac{2}{3})\). In the simplest 2-band model of graphene, each carbon atom carries one \(p_z\) orbital. We can add the orbitals by:

cell.add_orbital([1./3, 1./3], energy=0.0, label="pz")
cell.add_orbital([2./3, 2./3], energy=0.0, label="pz")

The fractional coordinate along c-axis is assumed to 0 if not specified. We can also make the orbitals to be located at \(z=0.5\times10\overset{\circ}{\mathrm {A}} = 5\overset{\circ}{\mathrm {A}}\) with:

cell.add_orbital([1./3, 1./3, 0.5], energy=0.0, label="pz")
cell.add_orbital([2./3, 2./3, 0.5], energy=0.0, label="pz")

The parameter energy gives the on-site energy of the orbital, which is assumed to be 0 eV if not specified. In absence of strain or external fields, the two orbitals have the same on-site energy. The parameter label is a tag to label the orbital.

In addition to fractional coordinates, the orbitals can also be added using Cartesian coordiantes with

# NOTE: a_1 should be parallel to x-axis
cell.add_orbital_cart([1.23, 0.71014083], unit=tb.ANG, energy=0.0, label="pz")
cell.add_orbital_cart([2.46, 1.42028166], unit=tb.ANG, energy=0.0, label="pz")

Here we use the parameter unit to specify the unit of the Cartesian coordinates. Note that PrimitiveCell class always uses fractional coordinates internally.

Add hopping terms

With the orbitals ready, we can add the hopping terms. The hopping terms for primitive cell with \(\gamma=60^\circ\) is shown in the figure below:

../../_images/graph_hop.png

Schematic plot of hopping terms of graphene. Primitive cells are indicated with dashed diamonds and numbered in blue text. Thick black arrows indicate the lattice vectors. Orbitals and hopping terms are shown as filled circles and red arrows, respectively.

From the figure we can see there are 6 hopping terms between \((0, 0)\) and neighbouring cells in the nearest approximation:

  • \((0, 0) \rightarrow (0, 0), i=0, j=1\)

  • \((0, 0) \rightarrow (0, 0), i=1, j=0\)

  • \((0, 0) \rightarrow (1, 0), i=1, j=0\)

  • \((0, 0) \rightarrow (-1, 0), i=0, j=1\)

  • \((0, 0) \rightarrow (0, 1), i=1, j=0\)

  • \((0, 0) \rightarrow (0, -1), i=0, j=1\)

Using the conjugate relation \(\langle i, 0 | \hat{H} | j, R\rangle = \langle j, 0 | \hat{H} | i, -R\rangle^*\) they can be reduced to:

  • \((0, 0) \rightarrow (0, 0), i=0, j=1\)

  • \((0, 0) \rightarrow (1, 0), i=1, j=0\)

  • \((0, 0) \rightarrow (0, 1), i=1, j=0\)

TBPLaS utilities the conjugate relation, so we need only to provide the reduced hopping terms. We do this by:

cell.add_hopping(rn=[0, 0], orb_i=0, orb_j=1, energy=-2.7)
cell.add_hopping(rn=[1, 0], orb_i=1, orb_j=0, energy=-2.7)
cell.add_hopping(rn=[0, 1], orb_i=1, orb_j=0, energy=-2.7)

The parameter rn specifies the index of neighbouring cell, while orb_i and orb_j give the indices of orbitals of the hopping term. energy is the hopping integral, which should be a complex number in general cases. The 3rd component of rn is assumed to be 0 if not provided. We can also specify it explicitly as:

cell.add_hopping(rn=[0, 0, 0], orb_i=0, orb_j=1, energy=-2.7)
cell.add_hopping(rn=[1, 0, 0], orb_i=1, orb_j=0, energy=-2.7)
cell.add_hopping(rn=[0, 1, 0], orb_i=1, orb_j=0, energy=-2.7)

Since monolayer graphene is two-dimensional, the 3rd component of rn is of cource 0. In general cases it should be non-zero.

Dump the model

Now we have successfully created the model. We can have a look at it by calling the plot method:

../../_images/graph_final.png

and print the orbitals and hopping terms by calling the print method. The output should look like:

Lattice vectors (nm):
    0.24600   0.00000   0.00000
    0.12300   0.21304   0.00000
    0.00000   0.00000   1.00000
Orbitals:
    0.00000   0.00000   0.00000 0.0
    0.33333   0.33333   0.00000 0.0
Hopping terms:
    (0, 0, 0) (0, 1) -2.7
    (1, 0, 0) (1, 0) -2.7
    (0, 1, 0) (1, 0) -2.7