I've only recently started using density matrices in my work but I am confused with the following code that I have whether I am getting the right matrix

def Hamiltonian(alpha,h):

    Sx = np.array([[0,1],[1,0]])
    Sy = np.array([[0,-1j],[1j,0]])
    Sz = np.array([[1,0],[0,-1]])
    I  = np.array([[1,0],[0,1]])

    H =  (alpha*np.kron(np.kron(Sx,Sx),I)) 
    H =+ (alpha*np.kron(np.kron(Sy,Sy),I)) 
    H =+ (alpha*np.kron(np.kron(I,Sx),Sx)) 
    H =+ (alpha*np.kron(np.kron(I,Sy),Sy)) 
    H =+ (h*np.kron(np.kron(I,Sz),I))

    return H

So the above gives me my Hamiltonian Function, where alpha is a real number and h is a magnetization paramteter acting on one of my qubits.

My ultimate goal with this is to plot |U| vs. h, where U = Tr(Hp). I know that the density matrix follows the form of the outer product of the wavefunction:

p = |ψ><ψ|

But my question is twofold:

1. How do I get the density matrix for the above Hamiltonian function? Has it got something to do with the eigenvalues and eigenvectors?

I have tried the following:

H = Hamiltonian(1,0.5)
print(H)

w, v = LA.eig(H)

print(w)
print(v)

ro = np.outer(w,w)
print(ro)

ro_square = ro**2
purity = np.trace(ro_square)
print(purity)

Comes out to be 0.5 which makes sense for my mixed state

U = (np.trace(H*ro))
print(U)

Alwyas U = 0 for all h...

But for every value of h I use I get U = 0 so something is going wrong, I feel like I'm not using the correct paramters for the density matrix in the outer product of the eigenvalues

2. How do I get the Hamiltonian in terms of an operator or is it already in this form?

Any help will be appreciated and if this is the wrong forum to put this question up let me know and I will put it up on the normal physics forum thanks.