Python Forum

My code is plotting eigenvalues of different H_T functions as E_c and E_v. For E_c i expect to have all E_c's for all Hamiltonians curved upward like for H_0. I can't find the mistake which cause it to go in wrong way for others. Belowe is my code

################################################################################
import math
import numpy as np
import matplotlib.pyplot as plt
import csv
n=150
alpha_up, alpha_down, beta_up, beta_down, kapa_up, kapa_down, Eta_up, Eta_down=-5.34, -5.71, -0.95, -0.52, -1.31, -1.23, 15.11, 17.10 #DFT
#alpha_up, alpha_down, beta_up, beta_down, kapa_up, kapa_down, Eta_up, Eta_down=-5.76, -6.2,-0.54, -0.03, -1.49, -1.4, 18.28, 20.93 #GW
#################################################################################################################################################################
f1, f2, a=1.4, 2.53/3.319, 3.319 #DFT
#f1, f2, a=2.26, 2.22/3.319, 3.319 #GW
##################################################################################################################################################################
Mp=0.25
K_p=0
Gamma=0.22
n=100
t=1
Qx=np.linspace(-Gamma/t,Mp/t,int(n))
Qy=0
###############################################################################################################################################################
#Pauli spin matricies
sx = np.array([[0, 1],[1, 0]], dtype=complex)
sy = np.array([[0, -1j],[1j, 0]], dtype=complex)
sz = np.array([[1, 0],[0, -1]], dtype=complex)
s1= np.array([[1,0],[0,1]], dtype=complex)
sg1=np.array([[1, 0], [0, 0]], dtype=complex)
sg2=np.array([[0, 0], [0, 1]], dtype=complex)
###############################################################################################################################################################
def H_T(k1,Alpha,Beta,Kapa,eta):
H_0=(0.5*f1*sz)+((f2*a)*((Qx[k1]*sx)+(Qy*sy)))
H_as=(Qx[k1]**2+Qy**2)*(Alpha*sg1+Beta*sg2)
H_3w=Kapa*(((Qx[k1]*sx)-(Qy*sy))**2)
H_cub=-(0.5*eta*(Qx[k1]**2+Qy**2))*((Qx[k1]*sx)+(Qy*sy))
H=H_0+H_as
if R==0: # H_0
H=H_0
elif R==1: #H_0+H_as
H=H_0+H_as
elif R==2: #H_0+H_as+H_3w
H=H_0+H_as+H_3w
elif R==3: #H_0+H_as+H_3w+H_cub
H=H_0+H_as+H_3w+H_cub
return H
#+H_as
#################################################################################################################################################################################################################################################################################################################
#Calculate eigenvalues of hamiltonian (E_c and E_v)
E_v=np.zeros(shape=(4,int(n)))
E_c=np.zeros(shape=(4,int(n)))
for R in range(4):
for j in range(n):
H1=H_T(j,alpha_up,beta_up,kapa_up,Eta_up)
e1=np.linalg.eigvals(H1).real
E_c[R,j]=e1[0]
H2=H_T(j,alpha_down,beta_down,kapa_down,Eta_down)
e2=np.linalg.eigvals(H2).real
E_v[R,j]=e2[1]
#Plot dispersion curve
fig1 = plt.figure()
plt.plot(Qx,E_c[0,:],'r',label=r'H$_0$', lw='3')
plt.plot(Qx,E_v[0,:], 'r--', lw='3')
plt.plot(Qx,E_c[1,:],'c',label=r'H$_0+$H$_{as}$', lw='3')
plt.plot(Qx,E_v[1,:], 'c--',lw='3')
plt.plot(Qx,E_c[2,:],'m',label=r'H$_0+$H$_{as}$+H$_{3w}$', lw='3')
plt.plot(Qx,E_v[2,:], 'm--',lw='3')
plt.plot(Qx,E_c[3,:],'k',label=r'H$_0$+H$_{as}$+H$_{3w}$+H$_{cub}$', lw='3')
plt.plot(Qx,E_v[3,:], 'k--', lw='3')
size=16
plt.legend(loc='best',fontsize=12)
plt.ylabel('Energy (eV)',fontsize=size)
plt.tick_params(axis = 'both', which = 'major', labelsize = 10)
fig1.savefig("Dispersion_MoSe2_Kormanyos15_DFT.pdf", bbox_inches='tight')
plt.show()
################################################################################

This post is so ugly.
reformat using BBcode (see buran's mod note for link to instructions) and repost.

if you call function H_T, it returns H,
so:

print('{}'.format(H_T(...))

will show results (note you supply arguments in place of ...