2
$\begingroup$

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.

enter image description here

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)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.