A first-principles calculation with SIESTA

The following is a very basic calculation of Z2 invariants using SIESTA for a generic system. Before running this calculation, an SCF run is required to create generate Hamiltonian and overlap matrices.

A complete example (including input files) for monolayer Bismuthen (111) can be found on GitHub.

import z2pack
import sisl
import numpy as np
import matplotlib.pyplot as plt

# Read geometry, Hamiltonian and overlap matrix from SIESTA output using sisl
sile = sisl.get_sile('in.fdf')
geom = sile.read_geometry()
H = sile.read_hamiltonian(geometry=geom)

# Create two helper functions to retrieve Hamiltonian and overlap for arbitrary k points
Hk = lambda k: H.Hk(k=k, dtype=np.complex64, format='array')
Sk = lambda k: H.Sk(k=k, dtype=np.complex64, format='array')

# Create array with the position of all orbitals
pos = geom.fxyz[list(map(geom.o2a, np.arange(H.no)))]
# In non-collinear/spin-orbit calculations each orbital is repeated twice once with up and once with down spin
pos = np.repeat(pos, 2, axis=0)

# Create system in z2pack
nEl = 0  # number of electrons (grep 'number of electrons' out.fdf)
system = z2pack.hm.System(

surface = lambda s, t: [t, s / 2, 0]

# Run the WCC calculations
settings = {
    # ...

result = z2pack.surface.run(system=system, surface=surface, **settings)