Berry flux over BZ

[bfocassio]

This is not much code related but it would serve as a nice example.

I'm trying to compute the Chern number for the Haldane model using the new feature that computes the berry flux and berry phase.

I'm setting up the model as:

graphene = sisl.geom.graphene(0.5774)
H = Hamiltonian(graphene,dtype=np.complex128)

delta=-0.2
t=-1.0
t2 =0.15*np.exp((1.j)*np.pi/2.)
t2c=t2.conjugate()

r = (0.1,  0.6)
t = (-delta , t )
H.construct([r, t])

H[1,1] = delta

H[0,0,(1,0)] = t2
H[1,1,(1,-1)] = t2
H[1,1,(0,1)] = t2
H[1,1,(1,0)] = t2c
H[0,0,(1,-1)] = t2c
H[0,0,(0,1)] = t2c

The bandstructure seems ok:

band = BandStructure(H, [[0.,0.],[2./3.,1./3.],[.5,0.5],[1./3.,2./3.], [0.,0.]], 400, [r'$\Gamma $',r'$K$', r'$M$', r'$K^\prime$', r'$\Gamma $'])

bs = band.asarray().eigh()

lk, kt, kl = band.lineark(True)
plt.xticks(kt, kl)
plt.xlim(0, lk[-1])
plt.ylim([-2.5, 2.5])
plt.ylabel(r'$E-\epsilon_F$ (eV)')
for bk in bs.T:
    plt.plot(lk, bk)

Now, for the berry flux, I'm doing:

nk = 24
kpts = np.linspace(-0.5,.5,nk,endpoint=False)

flux_lower_xy = np.zeros((len(kpts),len(kpts)))
for ki in range(len(kpts)):
    for kj in range(len(kpts)):
        origo = [kpts[ki], kpts[kj], 0.]
        flux_lower_xy[ki,kj] = H.eigenstate(k=origo).berry_flux()[0][0,1]

were the first index selects the lower band, and the last two the xy element.

And integrating it like:

simps([simps(flux_lower_xy[ki,:],kpts) for ki in range(len(kpts))],kpts)/(2*np.pi)

which should result in an integer. I'm obtaining something around 0.02 .

The model setup is ok ? And the berry flux part ?

Hoppings and berry flux/berry phase may be found here

Also, is there some way to open the boundary conditions (BC) of a periodic model (similar to the cut_piece of pythTB)?

Thanks in advance, any help is appreciated.

Version details 3.7.4 (default, Aug 13 2019, 20:35:49) [GCC 7.3.0] 0.9.8 13a327bd8e27d689f119bafdf38519bab7f6e0f6

[zerothi]

I think this is because my naming scheme for the functions are pretty bad... I am no expert in this, but I have tried to provide the necessary functionality.

These functions clearly need a revision to clarify their intent, thanks for bringing this to my attention.

[zerothi]

I have just implemented the Haldane model as Pythtb does it. There is however the difference that PythTB does not work in units the way sisl, does. So the two models are not exactly the same.

This may be checked by doing (see in the script about H.Hk(...)). One may add the same in PythTB (print(my_model._gen_ham([0.2, 0.2]))) but they are not the same, telling me that units are important here.

Note, that I haven't implemented the berry_flux and the name of the function was clearly wrong... :( Sorry about that.

If you have any ideas, please let me know.

import numpy as np
import sisl as si

import matplotlib.pyplot as plt

sc = si.SuperCell([[1, 0, 0], [0.5, 3 ** .5 / 2, 0], [0, 0, 10]], nsc=[3, 3, 1])
gr = si.Geometry(np.dot([[1/3, 1/3, 0], [2/3, 2/3, 0]], sc.cell), si.Atom(6, R=0.58), sc)

# System parameters
delta = 0.
t = -1.
t2 = 0.15 * np.exp(1j * np.pi / 2)
t2c = t2.conjugate()

H = si.Hamiltonian(gr, dtype=np.complex128)
H.construct([[0.1, .6], [0, t]])

# Create displaced on-site terms
H[0, 0] = -delta
H[1, 1] = delta

# Assing 2nd nearest neighbour couplings
# Contrary to PythTB sisl does not automatically impose
# symmetry. So you have to define all mirror couplings
# as well.
H[0, 0, (1, 0)] = t2
H[0, 0, (-1, 0)] = t2c
H[1, 1, (1,-1)] = t2
H[1, 1, (-1,1)] = t2c
H[1, 1, (0, 1)] = t2
H[1, 1, (0, -1)] = t2c
H[1, 1, (1, 0)] = t2c
H[1, 1, (-1, 0)] = t2
H[0, 0, (1, -1)] = t2c
H[0, 0, (-1, 1)] = t2
H[0, 0, (0, 1)] = t2c
H[0, 0, (0, -1)] = t2

# Pretty print-hoppings
#for r, c in H:
#    print("{r:2d} {c:2d} [ {i[0]:2d}, {i[1]:2d}]  == {v:.4f}".format(r=r, c=c % H.no, i=H.o2isc(c), v=H[r, c]))

# Compare the Hamiltonian with PythTB
#print(H.Hk([-0.2, -0.2, 0], format='array'))

special = [
    [0, 0, 0],
    [2/3, 1/3, 0],
    [1/2, 1/2, 0],
    [1/3, 2/3, 0],
    [0, 0, 0]]
label = [r'$\Gamma $',r'$K$', r'$M$', r'$K^\prime$', r'$\Gamma $']

band = si.BandStructure(H, special, 101, label)
lk, kt, kl = band.lineark(True)
plt.xticks(kt, kl)
plt.xlim(0, lk[-1])
plt.ylabel('$E-E_F$ [eV]')
plt.plot(lk, band.eigh())

# Now calculate the berry-phase
# Create a BZ model starting @ [-1/2, -1/2, 0]
Nx, Ny = 30, 30
mpx = si.MonkhorstPack(H, [Nx, 1, 1], trs=False)
mpy = si.MonkhorstPack(H, [1, Ny, 1], trs=False)
mpx.k[:, 0] = np.linspace(-0.5, 0.5, Nx, endpoint=False)
mpy.k[:, 1] = np.linspace(-0.5, 0.5, Ny, endpoint=False)

phi_a_1 = np.empty(Ny)
phi_b_1 = np.empty(Ny)
phi_c_1 = np.empty(Ny)

berry_phase = si.physics.electron.berry_phase
for ik, ky in enumerate(mpy):
    mpx.k[:, 1] = ky[1]
    # Calculate berry phase along line
    phi_a_1[ik] = berry_phase(mpx, sub=[0])
    phi_b_1[ik] = berry_phase(mpx, sub=[1])
    phi_c_1[ik] = berry_phase(mpx)

# I am not fully sure why this is needed?
phi_a_1 -= phi_a_1[0]
phi_b_1 -= phi_b_1[0]
phi_c_1 -= phi_c_1[0]

fig, ax = plt.subplots()
ky = mpy.k[:, 1]
ax.plot(ky, phi_a_1, 'ro')
ax.plot(ky, phi_b_1, 'go')
ax.plot(ky, phi_c_1, 'bo')
ax.set_title("Berry phase for lower (red), top (green), both bands (blue)")
ax.set_xlabel(r"$k_y$")
ax.set_ylabel(r"Berry phase along $k_x$")
ax.set_xlim(-0.5, 0.5)
ax.set_ylim(-7.,7.)
ax.yaxis.set_ticks([-2.*np.pi,-np.pi,0.,np.pi,2.*np.pi])
ax.set_yticklabels((r'$-2\pi$',r'$-\pi$',r'$0$',r'$\pi$', r'$2\pi$'))

plt.show()

[bfocassio]

This is awesome.

At the model setup, it is very important to note that we should impose the symmetry.

After setting up the model, if this is going to serve as an example, one could add a nice edge state visualization, such as:

pbcH = H.copy().tile(30,0)
obcH = pbcH.copy()
obcH.set_nsc(a=1)

ribbon_special = [[0.,0.,0.],[.0,.5,0.],[0.,0.,0.]]
ribbon_labels = [r'$\Gamma $', r'$Y$', r'$\Gamma $']

ribbon_band = si.BandStructure(pbcH, ribbon_special, 101, ribbon_labels)
lk, kt, kl = ribbon_band.lineark(True)

fig,ax = plt.subplots(1,2,sharex='all',sharey='all',figsize=(12,4))

ax[0].set_xticks(kt)
ax[0].set_xticklabels(kl)
ax[0].set_xlim(0, lk[-1])
ax[0].set_ylim([-2,2])
ax[0].set_ylabel('$E-E_F$ [eV]')

ax[0].plot(lk, ribbon_band.eigh())

ribbon_band.set_parent(obcH)
ax[1].plot(lk, ribbon_band.eigh())

ax[0].set_title('PBC')
ax[1].set_title('partial OBC')

plt.show()

At the berry phase part,

I would change it a bit to:

Nx, Ny = 31, 31
mpx = si.MonkhorstPack(H, [Nx, 1, 1], trs=False)
mpy = si.MonkhorstPack(H, [1, Ny, 1], trs=False)

After subtracting the first element (Note that the subtraction of the first element is needed to avoid the ambiguity of the berry phase definition).

If I'm not confusing things, the Chern number is simply given by the berry phase of occupied bands around the BZ contour, thus,

phi_a_1[-1]/(2*np.pi)

Is the berry phase routine working for SOC Hamiltonians ?

[zerothi]

At the model setup, it is very important to note that we should impose the symmetry.

Yes, I mention this one place in the documentation, but I should probably emphasize this.

After setting up the model, if this is going to serve as an example, one could add a nice edge state visualization, such as:

Thanks!!

pbcH = H.copy().tile(30,0)

you don't need to copy here, tile returns a new instance.

obcH = pbcH.copy()
obcH.set_nsc(a=1)

ribbon_special = [[0.,0.,0.],[.0,.5,0.],[0.,0.,0.]]
ribbon_labels = [r'$\Gamma $', r'$Y$', r'$\Gamma $']

ribbon_band = si.BandStructure(pbcH, ribbon_special, 101, ribbon_labels)
lk, kt, kl = ribbon_band.lineark(True)

fig,ax = plt.subplots(1,2,sharex='all',sharey='all',figsize=(12,4))

ax[0].set_xticks(kt)
ax[0].set_xticklabels(kl)
ax[0].set_xlim(0, lk[-1])
ax[0].set_ylim([-2,2])
ax[0].set_ylabel('$E-E_F$ [eV]')

ax[0].plot(lk, ribbon_band.eigh())

ribbon_band.set_parent(obcH)
ax[1].plot(lk, ribbon_band.eigh())

ax[0].set_title('PBC')
ax[1].set_title('partial OBC')

plt.show()

At the berry phase part,

I would change it a bit to:

Nx, Ny = 31, 31
mpx = si.MonkhorstPack(H, [Nx, 1, 1], trs=False)
mpy = si.MonkhorstPack(H, [1, Ny, 1], trs=False)

It does not matter, a MP grid never uses the same point twice, so for Nx = 4, you would get [-1/4, 0, 1/4, 1/2]

After subtracting the first element (Note that the subtraction of the first element is needed to avoid the ambiguity of the berry phase definition).

If I'm not confusing things, the Chern number is simply given by the berry phase of occupied bands around the BZ contour, thus,

phi_a_1[-1]/(2*np.pi)

No this won't work (I guess), you would have to use the first and last point to get the accumulated Berry-phase over the path (this is how I currently understand it, but I could be wrong!)

Is the berry phase routine working for SOC Hamiltonians ?

I think so, the Berry phase only uses the eigenstates. However, the berry_curvature will not work since it requires the velocity matrix.

[bfocassio]

Thank you for your patience and help.

No this won't work (I guess), you would have to use the first and last point to get the accumulated Berry-phase over the path (this is how I currently understand it, but I could be wrong!)

You're right, the best way would be (phi_a_1[0] - phi_a_1[-1])/(2*np.pi)

Now I'm going to explore the Kane-Mele model

[zerothi]

Hmm, shouldn't it be phi_a_1[-1] - phi_a_1[0] since you want the accumulated phase? if the last is 2pi above the first, it would give you Chern == 1

[bfocassio]

Yes, that's what I meant, thanks again

[bfocassio]

I set up the Kane and Mele model (as in here) and I'm trying to calculate the Berry phase for the two occupied bands, together and separately.

The model is:

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

# ### Setting up the model

sc = si.SuperCell([[1., 0, 0], [0.5, 3 ** .5 / 2, 0], [0, 0, 30]], nsc=[3, 3, 1])
gr = si.Geometry(np.dot([[1/3, 1/3, 0], [2/3, 2/3, 0]], sc.cell), si.Atom(6, R=0.58), sc)

# ### Useful definitions
#[Re(upup) Re(dndn) Re(updn) Im(updn) Im(upup) Im(dndn) Re(dnup) Im(dnup)]

eye = np.array([1.,1.,0.,0.,0.,0.,0.,0.])
sigma_z = np.array([1,-1,0,0,0,0,0,0])
sigma_x = np.array([0,0,1,0,0,0,1,0])
sigma_y = np.array([0,0,0,-1,0,0,0,1])

imag_sigma_z = np.array([0,0,0,0,1.,-1.,0,0])
imag_sigma_x = np.array([0,0,0,1,0,0,0,1])
imag_sigma_y = np.array([0,0,1,0,0,0,-1,0])

# System parameters

t = 1.
esite = .1*t
spin_orb = .06*t
rashba = .05*t

H = si.Hamiltonian(gr,dim=8)

# spin-independent first-neighbor hoppings

H.construct([[0.1, .6], [0, t*eye]])

# set on-site energies

H[0, 0] = esite*eye
H[1, 1] = -1*esite*eye

# Assing 2nd nearest neighbour couplings

# Define all mirror couplings as well.

# second-neighbour spin-orbit hoppings (s_z)

H[0, 0, (0,1)] = -1*spin_orb*imag_sigma_z
H[0, 0, (0,-1)] = 1*spin_orb*imag_sigma_z

H[0, 0, (1,0)] = 1*spin_orb*imag_sigma_z
H[0, 0, (-1,0)] = -1*spin_orb*imag_sigma_z

H[0, 0, (1,-1)] = -1*spin_orb*imag_sigma_z
H[0, 0, (-1,1)] = 1*spin_orb*imag_sigma_z

H[1, 1, (0,1)] = 1*spin_orb*imag_sigma_z
H[1, 1, (0,-1)] = -1*spin_orb*imag_sigma_z

H[1, 1, (1,0)] = -1*spin_orb*imag_sigma_z
H[1, 1, (-1,0)] = 1*spin_orb*imag_sigma_z

H[1, 1, (1,-1)] = 1*spin_orb*imag_sigma_z
H[1, 1, (-1,1)] = -1*spin_orb*imag_sigma_z

# Rashba first-neighbor hoppings: (s_x)(dy)-(s_y)(d_x)

r3h =np.sqrt(3.0)/2.0

H[0,1, (0,0)] +=  1*rashba*(0.5*imag_sigma_x-r3h*imag_sigma_y)
H[1,0, (0,0)] += -1*rashba*(0.5*imag_sigma_x-r3h*imag_sigma_y)

H[0,1, (0,-1)] += 1*rashba*(-1.0*imag_sigma_x)
H[0,1, (0,1)] += -1*rashba*(-1.0*imag_sigma_x)

H[0,1, (-1,0)] += 1.*rashba*( 0.5*imag_sigma_x+r3h*imag_sigma_y)
H[0,1, (1,0)] += -1.*rashba*( 0.5*imag_sigma_x+r3h*imag_sigma_y)

# print model at gamma

print(H.Hk([0., 0., 0], format='array'))

# bulk band structure

special = [
    [0, 0, 0],
    [2/3, 1/3, 0],
    [1/3, 2/3, 0],
    [0, 0, 0]]
label = [r'$\Gamma $',r'$K$', r'$K^\prime$', r'$\Gamma $']

band = si.BandStructure(H, special, 301, label)
lk, kt, kl = band.lineark(True)
plt.xticks(kt, kl)
plt.xlim(0, lk[-1])
plt.ylabel('$E-E_F$ [eV]')
plt.plot(lk, band.eigh())
plt.show()

# Ribbon band structure for partially finite geometry

pbcH = H.tile(30,0)

obcH = pbcH.copy()
obcH.set_nsc(a=1)

ribbon_special = [[0.,0.,0.],[.0,.5,0.],[0.,0.,0.]]
ribbon_labels = [r'$\Gamma $', r'$Y$', r'$\Gamma $']

ribbon_band = si.BandStructure(pbcH, ribbon_special, 501, ribbon_labels)
lk, kt, kl = ribbon_band.lineark(True)

fig,ax = plt.subplots(1,2,sharex='all',sharey='all',figsize=(12,4))

ax[0].set_xticks(kt)
ax[0].set_xticklabels(kl)
ax[0].set_xlim(0, lk[-1])
ax[0].set_ylim([-1.75,1.75])
ax[0].set_ylabel('$E-E_F$ [eV]')

ax[0].plot(lk, ribbon_band.eigh())

ribbon_band.set_parent(obcH)
ax[1].plot(lk, ribbon_band.eigh())

ax[0].set_title('PBC')
ax[1].set_title('partial OBC')

plt.show()

For the Berry phase I do:

# Now calculate the berry-phase

# Create a BZ model starting @ [-1/2, -1/2, 0]

Nx, Ny = 30,30
mpx = si.MonkhorstPack(H, [Nx, 1, 1], trs=False)
mpy = si.MonkhorstPack(H, [1, Ny, 1], trs=False)
mpx.k[:, 0] = np.linspace(-0.5, 0.5, Nx, endpoint=False)
mpy.k[:, 1] = np.linspace(-0.5, 0.5, Ny, endpoint=False)

phi_occ_0 = np.empty(Ny)
phi_occ_1 = np.empty(Ny)
phi_occ = np.empty(Ny)

# Calculate berry phase along line

berry_phase = si.physics.electron.berry_phase
for ik, ky in enumerate(mpy):
    mpx.k[:, 1] = ky[1]

    phi_occ_0[ik] = berry_phase(mpx, sub=[0])
    phi_occ_1[ik] = berry_phase(mpx, sub=[1])
    phi_occ[ik] = berry_phase(mpx, sub=[0,1])

This returns me the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-39-4fd182f2dfe7> in <module>
      3     mpx.k[:, 1] = ky[1]
      4     # Calculate berry phase along line
----> 5     phi_occ_0[ik] = berry_phase(mpx, sub=[0])
      6     phi_occ_1[ik] = berry_phase(mpx, sub=[1])
      7     phi_occ[ik] = berry_phase(mpx, sub=[0,1])

~/anaconda3/lib/python3.7/site-packages/sisl/physics/electron.py in berry_phase(contour, sub, eigvals, closed, method)
   1188 
   1189     # Do the actual calculation of the final matrix
-> 1190     d = _berry(contour.asyield().eigenstate())
   1191 
   1192     # Correct return values

~/anaconda3/lib/python3.7/site-packages/sisl/physics/electron.py in _berry(eigenstates)
   1169         def _berry(eigenstates):
   1170             first = next(eigenstates).sub(sub)
-> 1171             first.change_gauge('r')
   1172             prev = first
   1173             prd = 1

~/anaconda3/lib/python3.7/site-packages/sisl/physics/electron.py in change_gauge(self, gauge)
   1796 
   1797         if gauge == 'r':
-> 1798             self.state *= np.exp(1j * phase).reshape(1, -1)
   1799         elif gauge == 'R':
   1800             self.state *= np.exp(-1j * phase).reshape(1, -1)

ValueError: operands could not be broadcast together with shapes (1,4) (1,2) (1,4) 

Any ideas of what I'm doing wrong ?

Thanks in advance

[zerothi]

Thanks for discovering this... It is definitely a bug!

[zerothi]

@bfocassio could you try the latest commit of sisl?

[bfocassio]

Yes, It's working

Thank you

[zerothi]

@bfocassio I am going to re-open.

Then I will remember to make a small tutorial about (if you don't mind me using your supplied codes?)

:)

[bfocassio]

I don't mind at all. And thank you again