How Does the Metropolis Algorithm Calculate Specific Heat in a 2D XY Model?

[Firben]

Homework Statement

In a XY - model in two dimensional square lattice. I want to calculate specific heat capacity, <E> and <E^2>. By using the formula:

H = -Jsum(<i,j>)Cos(tetha(i)-tetha(j))

The Metropolis algorithm reads:

Simulations for Statistical and Thermal Physics

Homework Equations

H = -Jsum(<i,j>)Cos(tetha(i)-tetha(j))

The Attempt at a Solution

Here is my code for the Hamiltonian:
http://s716.photobucket.com/user/Pitoraq/media/MC_zpsqqzwctgg.png.html?sort=3&o=0

What is wrong with my code ?

 

[Firben]

Updated my code:

clear
clc
J = 1;
for N = 1:10000

for M = 1:2
tetha1 = (-1 + (1+1)*rand(4,1))*2*pi;
tetha2 = (-1 + (1+1)*rand(4,1))*2*pi;
for i=1:1:4

Hamton(i,M) = J*(cos((tetha1(i)-tetha2(i))));

end

Hamtonsum = sum(Hamton);

end
DeltaE(N) = Hamtonsum(1) - Hamtonsum(2);
end

Still clueless what to do

 

[magoo]

Is pi defined anywhere?

 

[Firben]

Pi is build in Matlab, no need to define it

 

[magoo]

OK, I've done some C++ but not MatLab.

How about Hamtonsum, Hamtonsum(1) and Hamtonsum(2) ? Should Hamtonsum have a dimension?

 

[Firben]

That is not the problem now. It is Hamton(i,M) = J*(cos((tetha1(i)-tetha2(i)))); .
Why should i take the hamiltonian 8 times ?. If tetha = 2*pi*[-1,1] (random numbers).
H = -(cos(tehta(i) - tetha(j) + cos(tetha(i) - etha(j+1)+..)
How do i choose the value tetha(i) ?

 

[Firben]

My new code is :

clear
clc
Nsq = 100*100;
N = sqrt(Nsq);
T = linspace(0,2,100);J = 1;

randIndx = randi([1,Nsq],Nsq,1);
randR = rand(Nsq,1);
for j = 1:100
spins=2*pi*rand(N);

for i=1:Nsq
index = randIndx(i);
H = 0;
H = H + cos(spins(mod(index,Nsq) + 1) - spins(mod(index - 1 + N,Nsq) + 1));
H = H + cos(spins(mod(index - 2 + Nsq,Nsq) + 1) - spins(mod(index - 1 + Nsq,Nsq) + 1));
H = -J*H;
%Hamiltonian(j) = -H; Hi = (spins(index)*H)/2;
if(Hi <= 0 || randR(i) < (spins(index)*H)/2 )
spins(index) = -spins(index);
end

E(j) = (1/j)*H;

%plot(T,E,'r');
%plot(T,H);
end

%plot(T,Cv,'b');
end
%Esquared(j) = (1/j).*H.^2
%Cv = (1./T.^2).*((E.^2)-Esquared);

plot(T,E);

what is wrong ?