# Homework Help: Metropolis algorithm in MC

1. Feb 21, 2017

### Firben

1. The problem statement, all variables and given/known data
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))

Simulations for Statistical and Thermal Physics

2. Relevant equations
H = -Jsum(<i,j>)Cos(tetha(i)-tetha(j))

3. 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 ?

2. Feb 23, 2017

### 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

3. Feb 23, 2017

### magoo

Is pi defined anywhere?

4. Feb 23, 2017

### Firben

Pi is build in Matlab, no need to define it

5. Feb 23, 2017

### magoo

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

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

6. Feb 27, 2017

### 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) ?

7. Mar 15, 2017

### 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 ?