# Heisenberg Model python : Specific heat capacity for spin 2

I have the correct plot for specific heat capacity when I am using the formula which is $$C_V$$ = differentiation of entropy with respect to temperature. However, When I try to calculate $$C_V$$, by using formula

$$C_V = \frac{\langle E^2\rangle - \langle E\rangle ^2}{T^2}\, ,$$

I am getting a strange result. I am expecting the second figure just like the first one. Specific heat doesn't need to go negative however it has many negative values in the second figure.

Following is the code, any help would be immensely welcome.

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(789645)
# CONSTANTS INSTANTIATIONS          ----------------------------------------------------
J =1.0  ;Jeff = 2*J;D = 0.*Jeff;h = 0;Spin = 2
n = int(2*Spin + 1)  #size of the matrix
k = Spin*(Spin+1)   # total value spin can take
T1 = 0.1; T2 = 5;S = np.arange(-Spin,Spin+1e-5)   #possible values of S
mm=30;T = np.linspace(0.1,6.,mm);beta = [1./a for a in T]
#---------------------------------------------------
Su = np.zeros((n,n));Sd = np.zeros((n,n));Sz = np.zeros((n,n))

np.fill_diagonal(Sz,S)   # this gives Sz which is diagonal matrix. we
#Assume   Sz to be the same for both the lattice A and B
for i in range(n-1):
Su[i,i+1] = np.sqrt(k - (S[i]+1)*S[i])
Sd[i+1,i] = np.sqrt(k - (S[i]+1)*S[i])

Sx=0.5*(Su+Sd);Sy=(-0.5j)*(Su-Sd)

def ham(mA,mB,site):
if site == "A":
Ham =-D*Sz*Sz-h*Sz+Jeff*(mB['x']*Sx+mB['z']*Sz)- \
Jeff/2.0*(mA['x']*mB['x'] + mA['z']*mB['z'])*np.eye(n)
else:
Ham =-D*Sz*Sz-h*Sz+Jeff*(mA['x']*Sx+mA['z']*Sz)-\
Jeff/2.0* (mA['x']*mB['x'] + mA['z']*mB['z'])*np.eye(n)
return np.linalg.eigh(Ham)

def expectation(b,en,ev,spin):
Z=np.sum(np.exp(-b*en))
ave = 0.
for i in range(n):
ave += np.exp(-b*en[i])*(np.matmul(np.asmatrix(np.reshape \
(ev[:,i],(n,1))).getH(),np.matmul(np.asmatrix(spin), \
np.asmatrix(np.reshape(ev[:,i],(n,1))))))
return ave[0,0]/Z

def FreeEnergy(b,en):
Z=np.sum(np.exp(-b*en))
FE=-1/b*np.log(Z)
return (FE,Z)

nn = 1   # samplin number
ff = []  # Free energy save here
mA_list = []  # expectation of <S>A
mB_list = []  # expectation of <S>B
MxA = [];MzA = [];MxB = [];MzB = [];Mx = [];Mz = [];ent = [];Cv = []
mA = {'x':np.random.rand(),'z':np.random.rand()}
mB = {'x':np.random.rand(),'z':np.random.rand()}

for bt in beta:
change=1
itr=0
while change > 1e-5:

itr += 1
mA_new={}
mB_new={}

enA,evA = ham(mA,mB,"A")
mA_new['x']=expectation(bt,enA,evA,Sx)
mA_new['z']=expectation(bt,enA,evA,Sz)

enB,evB = ham(mA,mB,"B")
mB_new['x']=expectation(bt,enB,evB,Sx)
mB_new['z']=expectation(bt,enB,evB,Sz)
change=abs(mA_new['x']-mA['x'])+\
abs(mA_new['z']-mA['z'])+\
abs(mB_new['x']-mB['x'])+\
abs(mB_new['z']-mB['z'])
p=0.2
mA['x']=p*mA['x']+(1-p)*mA_new['x']
mA['z']=p*mA['z']+(1-p)*mA_new['z']
mB['x']=p*mB['x']+(1-p)*mB_new['x']
mB['z']=p*mB['z']+(1-p)*mB_new['z']

mA_list.append((mA['x'],mA['z']))
mB_list.append((mB['x'],mB['z']))
(FreeA,ZA)=FreeEnergy(bt,enA)
(FreeB,ZB)=FreeEnergy(bt,enB)

# These are the average energy  <E>
UA = np.sum(np.dot(enA,np.exp(-bt*enA)))/ZA
UB = np.sum(np.dot(enB,np.exp(-bt*enB)))/ZB
UT = UA + UB

ent.append((UT-FreeA-FreeB)*bt)   # This gives the entropy, F = U - TS

#  These are <E^2>, I think the problem lies here.
UA2 = np.sum(np.dot(enA**2,np.exp(-bt*enA)))/ZA
UB2 = np.sum(np.dot(enB**2,np.exp(-bt*enB)))/ZB
UT2 = UA2 + UB2
Cv.append(((UT2 - UT**2))*bt**2)

sp_heat = np. diff(ent, n =1)
deltaT = (T2-T1)/mm
sp_heat[:] = T[:mm-1] * [a/deltaT for a in sp_heat]

# The first plot is correct. Which is obtained from the derivative of
# Entropy, where the second plot is incorrect
# and should appear as first is obtained from the fluctuation of energy
#<E^2> -<E>^2)*bt**2

#plt.plot(T[:mm-1], sp_heat)

plt.plot(T, Cv)