sisl
sisl copied to clipboard
zerothi
Berry flux over BZ
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
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.
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()
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 ?
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.
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
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
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
@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?)
:)