from __future__ import division, print_function import numpy as np import matplotlib import kwant from ipywidgets import interact #from IPython.html.widgets import interact #from ipywidgets import StaticInteract, RangeWidget, DropDownWidget #from IPython.display import display_html from matplotlib import pyplot as plt from mpl_toolkits.mplot3d import Axes3D import pfapack as pf #from edx_components import * #%run ../code/init_mooc_nb.ipy from matplotlib import cm #from ipywidgets import RadioWidget import scipy class SimpleNamespace(object): def __init__(self, **kwargs): self.__dict__.update(kwargs) def haldane_ribbon(w, boundary='zigzag'): def ribbon_shape_zigzag(pos): fac = np.sqrt(3)/2 (x, y) = pos return (-0.5 / np.sqrt(3) - 0.1 <= y < fac * w + 0.01 ) def ribbon_shape_armchair(pos): fac = np.sqrt(3)/2 (x, y) = pos return (-1 <= x < w ) def onsite(site, par): if site.family == a: return par.m else: return -par.m def nn_hopping(target, source, par): return par.t def nnn_hopping(target, source, par): return 1j*par.t_nnn graphene = kwant.lattice.honeycomb() a, b = graphene.sublattices nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a)) nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b)) if boundary == 'zigzag': sym = kwant.TranslationalSymmetry((1, 0)) else: sym = kwant.TranslationalSymmetry((0, np.sqrt(3))) sys = kwant.Builder(sym) if boundary == 'zigzag': sys[graphene.shape(ribbon_shape_zigzag, (0, 0))] = onsite else: sys[graphene.shape(ribbon_shape_armchair, (0, 0))] = onsite sys[graphene.neighbors()] = nn_hopping sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = nnn_hopping sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = nnn_hopping return sys.finalized() #returns a kwant system, that resembles the momentum subblock of the Hamiltonian def Haldane_infinite(): lat = kwant.lattice.chain() sys = kwant.Builder() def onsite0(site, par): phases_nnn = np.array([np.sqrt(3) * par.ky, 1.5 * par.kx - np.sqrt(3) / 2 * par.ky, -1.5 * par.kx - np.sqrt(3) / 2 * par.ky]) temp = np.sum(np.exp(1j * phases_nnn + 1j * par.phi)) temp = temp + temp.conjugate() return par.m + par.t_nnn * temp def onsite1(site, par): phases_nnn = np.array([np.sqrt(3) * par.ky, 1.5 * par.kx - np.sqrt(3) / 2 * par.ky, -1.5 * par.kx - np.sqrt(3) / 2 * par.ky]) temp = np.sum(np.exp(-1j * phases_nnn + 1j * par.phi)) temp = temp + temp.conjugate() return - par.m + par.t_nnn * temp def hopping (site1, site2, par): phases = np.array([par.kx, (np.sqrt(3) * par.ky - par.kx) / 2, -(np.sqrt(3) * par.ky + par.kx)/2]) return par.t * np.sum(np.exp(1j * phases)) sys[lat(0)] = onsite0 sys[lat(1)] = onsite1 sys[kwant.HoppingKind((1,), lat)] = hopping return sys.finalized() def Qi_Wu_Zhang_infinite(): sigma_x = np.array([[0, 1], [1, 0]]) sigma_y = np.array([[0, -1j], [1j, 0]]) sigma_z = np.array([[1, 0], [0, -1]]) def onsite(site, par): return (par.delta * np.sin(par.kx) * sigma_y + 2 * par.gamma * np.sin(par.ky) * sigma_x + (2 * par.gamma * np.cos(par.ky) + par.mu + 2 * par.t * np.cos(par.kx)) * sigma_z) lat = kwant.lattice.chain() sys = kwant.Builder() sys[lat(0)] = onsite return sys.finalized() # dispersion functions def plot_2D(X,Y,Z, ax_in=None): if ax_in==None: fig = plt.figure(figsize=(7,5)) ax = fig.add_subplot(111, projection='3d') else: ax = ax_in vmin = np.array(Z).min() vmax = np.array(Z).max() if len(np.shape(Z)) > 2: for z in Z: ax.plot_surface(X, Y, z, rstride=1, cstride=1, cmap=cm.RdBu_r, linewidth=0.1, vmin=vmin, vmax=vmax) else: ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.RdBu_r, linewidth=0.1, vmin=vmin, vmax=vmax) if ax_in==None: return fig, ax else: return ax def evaluate_on_grid(X, Y, func): """ X, Y should be in np.meshgrid form. It's enough for func to work on floats. """ data = [] for xx, yy in zip(X, Y): row = [] for i,j in zip(xx, yy): row.append(func(i,j)) data.append(row) data = np.array(data) data = [np.array(data[:,:,i]) for i in range(np.shape(data)[2])] return data def diagonalize(sys, par): mat = sys.hamiltonian_submatrix(args=[par]) ev, evec = scipy.linalg.eigh(mat) # Automatically sorted return ev, evec def dispersion_func(sys, par): def func(kx, ky): par.kx = kx par.ky = ky return diagonalize(sys, par)[0] return func def eigen_func(sys, par): def func(kx, ky): par.kx = kx par.ky = ky return diagonalize(sys, par)[1] return func def plot_haldane_dispersion(t_nnn, m=0): par = SimpleNamespace(t=1.0, t_nnn=t_nnn, m=m, phi=np.pi/2, kx=0.0, ky=0.0) sys = Haldane_infinite() K = np.linspace(-np.pi, np.pi, 150) mesh = np.meshgrid(K, K) energies = evaluate_on_grid(*mesh, func=dispersion_func(sys, par)) fig, ax = plot_2D(*mesh, Z=energies) ax.set_xlabel('$k_x$') ax.set_xticks([-np.pi, 0.0, np.pi]) ax.set_xticklabels(['$-\pi$', '$0$', '$\pi$']) ax.set_ylabel('$k_y$') ax.set_yticks([-np.pi, 0.0, np.pi]) ax.set_yticklabels(['$-\pi$', '$0$', '$\pi$']) ax.set_zlabel('$E$') ax.set_zticks([-3, 0.0, 3]) ax.set_zticklabels(['$3$', '$0$', '$3$']) ax.view_init(8,20) ax.set_zlim3d(-4,4) return fig def berry_curvature(eigen, num_points=51, num_filled_bands=1): K = np.linspace(-np.pi, np.pi, num_points, endpoint=False) vectors = np.array([[eigen(kx, ky)[:, :num_filled_bands] for kx in K] for ky in K]) vectors_x = np.roll(vectors, 1, 0) vectors_xy = np.roll(vectors_x, 1, 1) vectors_y = np.roll(vectors, 1, 1) shifted_vecs = [vectors, vectors_x, vectors_xy, vectors_y] v_shape = vectors.shape shifted_vecs = [i.reshape(-1, v_shape[-2], v_shape[-1]) for i in shifted_vecs] dets = np.ones(len(shifted_vecs[0]), dtype=complex) for vec, shifted in zip(shifted_vecs, np.roll(shifted_vecs, 1, 0)): dets *= [np.linalg.det(a.T.conj().dot(b)) for a, b in zip(vec, shifted)] bc = np.angle(dets).reshape(int(np.sqrt(len(dets))), -1) bc = (bc + np.pi/2) % (np.pi) - np.pi/2 return bc plot_haldane_dispersion(0.0,0)