In [5]:

    

from sympy import symbols, pi, I,mpmath,cos, sin
from sympy.matrices import Matrix
from sympy.interactive.printing import init_printing
from sympy.physics.matrices import msigma
from sympy.functions import conjugate
from customOperators import*
import numpy as np
from numpy import linalg as LA
import math as math
from scipy.linalg import eig
#import scipy.sparse as sparse
from scipy.linalg import eigvals


np.set_printoptions(suppress=True)


init_printing()

pi=pi.evalf(40)
A,B,C,D =symbols('A:D')
R_0 =symbols('R_0')
a = symbols('a')
M=symbols('M')
k=symbols('k')
tA= A/(a*2)
tB=B/(a**2)
tD=D/(a**2)


I_2 =  Matrix(2, 2, lambda i,j: int(i==j))


def ctranspose(matrix):
    return matrix.transpose().applyfunc(lambda i:conjugate(i))

class Hamiltonian(object):
    
    def __init__(self,nY):
        self.H_0 = (C-float(4)*tD)*(I_2|kron|I_2) +(M-float(4)*tB)*(I_2|kron|msigma(3)) 
        #I*A/a*((msigma(3)|kron|msigma(1))-(I_2|kron|msigma(2)))
        #I*R_0/a*(TensorProduct(msigma(2),(msigma(3)+I_2)/2)-TensorProduct(msigma(1),(msigma(3)+I_2)/2))
        
        self.V_0=  +(tD)*(I_2|kron|I_2)+(tB)*(I_2|kron|msigma(3))-\
        I*tA*(I_2|kron|msigma(2))
        #I*R_0/a*TensorProduct(msigma(1),(msigma(3)+I_2)/2)
        
        self.W_0d =  +(tD)*(I_2|kron|I_2)+(tB)*(I_2|kron|msigma(3))
        self.W_0off = -I*tA*(msigma(3)|kron|msigma(1))
        #I*A/(2*a)*(msigma(3)|kron|msigma(1))
        #I*R_0/a*TensorProduct(msigma(2),(msigma(3)+I_2)/2) 
        
        self.IdDown = Matrix(nY,nY , lambda i,j :int(i-j==1))
        self.IdUp = Matrix(nY,nY , lambda i,j :int(i-j==-1))
        self.Id= Matrix(nY,nY , lambda i,j :int(i==j))
    

    def hamiltonianNumeric(self, valA,valB,valC,valD,valR_0,vala,valM,valk):
        matrix= self.expandHilbertSpace().evalf(subs={A:valA,B:valB,C:valC,D:valD,R_0:valR_0,a:vala,M:valM,k:valk})
        return  np.array(matrix.tolist()).astype(np.float64)
    
    def V(self, valA,valB,valC,valD,valR_0,vala,valM):
        matrix =  (self.Id|kron|self.H_0) +\
               (self.IdUp|kron|self.V_0)+(self.IdDown|kron|ctranspose(self.V_0))
        numeric = matrix.evalf(subs={A:valA,B:valB,C:valC,D:valD,R_0:valR_0,a:vala,M:valM})
        return {
            'numeric':np.array(numeric.tolist(),dtype=float),
            'object':matrix
        }
    
    def Wcos(self, valA,valB,valC,valD,valR_0,vala,valM):
        matrix=  2*((self.Id|kron|self.W_0d))
        numeric = matrix.evalf(subs={A:valA,B:valB,C:valC,D:valD,R_0:valR_0,a:vala,M:valM})
        return {
            'numeric':np.array(numeric.tolist(),dtype=float),
            'object':matrix
        }    
    
    def Wsin(self, valA,valB,valC,valD,valR_0,vala,valM):
        matrix=  2*I*(self.Id|kron|self.W_0off)
        numeric = matrix.evalf(subs={A:valA,B:valB,C:valC,D:valD,R_0:valR_0,a:vala,M:valM})
        return {
            'numeric':np.array(numeric.tolist(),dtype=float),
            'object':matrix
        }
    
    def eigenvals(self,k,matrixV,matrixWcos,matrixWsin):

        k=k*pi
        seno= math.sin(k)
        cosseno= math.cos(k)
        w= eigvals(np.array(matrixV+matrixWcos*cosseno+matrixWsin*seno))
        w=np.asarray(w,np.float64)
        idx = w.argsort()[::-1]   
        w = w[idx]
        return w


    

In [ ]: