# BHZ model¶

This example shows the calculation of the $$\mathbb{Z}_2$$ invariant for the periodic BHZ model.

import logging

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

logging.getLogger('z2pack').setLevel(logging.WARNING)

pauli_y = np.array([[0, 1], [1, 0]], dtype=complex)
pauli_x = np.array([[0, -1j], [1j, 0]])
pauli_z = np.diag([1, -1])
pauli_vector = [pauli_x, pauli_y, pauli_z]

def bhz(A, B, C, D, M):
def h(k):
kx, ky = 2 * np.pi * np.array(k)
d = [
A * np.sin(kx), -A * np.sin(ky),
-2 * B * (2 - (M / (2 * B)) - np.cos(kx) - np.cos(ky))
]
H = sum(di * pi for di, pi in zip(d, pauli_vector))
epsilon = C - 2 * D * (2 - np.cos(kx) - np.cos(ky))
H += epsilon
return H

def Hamiltonian(k):
return np.vstack([
np.hstack([h(k), np.zeros((2, 2))]),
np.hstack([np.zeros((2, 2)),
h(-np.array(k)).conjugate()])
])

return Hamiltonian

def run(**kwargs):
system = z2pack.hm.System(bhz(**kwargs), dim=2)
res = z2pack.surface.run(system=system, surface=lambda s, t: [s / 2, t])
z2pack.plot.wcc(res)
plt.savefig('wcc_plot.pdf', bbox_inches='tight')
print('Z2 invariant:', z2pack.invariant.z2(res))

if __name__ == '__main__':
run(A=0.5, B=1., C=0., D=0., M=1.)