Take the 2-minute tour ×

i have a piece of metropolis algorithm:

mB=5.79*10^(-9); %Bohr magnetone in eV*G^-1
kB=0.86*10^(-4); %Boltzmann in eV*K^-1
%system parameters
L=60; %side square grid
L2=L*L; % total number grid position
Tstep=5; %step in temperature change (K)
Maxstep=10; %max number of steps 
nmcs=5; % cycle numberof Metropolis algorithm

magnet=NaN(1,Maxstep);%store magnetization in "monte carlo images" of sample

%Creation initial point arrangement of magnetic spins
%Outer parameters
H=100000; %Gauss
T=20; % Kelvin

%Energy alteration in spin-reverse
de =@ (i,j) (2*mB*H).*mlat(i,j);
%Metropolis probability
pmetro=@ (i,j) exp(-de(i,j)./(kB*T));

%Creation and display of initial lattice


% Alteration of system with time
for ii=1:Maxstep
    for imc=1:nmcs
        for i=1:L
            for j=1:L
                 if pmetro(i,j)>=1

                elseif rand<pmetro(i,j)

   % shading interp;

error=std(magnet) ./sqrt(numel(magnet))
fprintf('Temperature = %d K',T)

axis([0 10 0 0.5])
grid on
xlabel('i (Configuration) ')

Now,the problem is in figure(13).The values it gives me are around zero (0.05,0.02..).It supposes to give me values around 0.3.. Generally,the graph its ok,It gives me the right "shape"(it has points) but as i said around zero. I really don't know how to put this post in order to be understood.Maybe i have some mistake in the "magnet"matrix ,i don't know. Anyway,i don't demand from anybody to check it thoroughly ,i am just asking if with a quick look anyone can help.

ΕDIT--->> Also,sometimes when i run the program ,it gives me :

Undefined function or method 'mlat' for input arguments of type 'double'.

Error in ==> @(i,j)(2*mB*H).*mlat(i,j)

Error in ==> @(i,j)exp(-de(i,j)./(kB*T))

Error in ==> metropolis at 39 if pmetro(i,j)>=1

EDIT--->>> I found the "mistake" .In my code in the loops where i have the function "pmetro" i replaced it with the "exp(-(2*mB*H).*mlat(i,j)./(kB*T))" and the program worked just fine!!! Why it didn't work with calling the "pmetro"??How can i overcome this?Is there a problem with function handles in loops?


share|improve this question
-1: This community of volunteers offers advice and suggestions. Being asked to debug this code is like asking us to do your work for you. Please ask a precise question and re-submit, we'll happily try to help! –  André Caron Feb 22 '11 at 16:14
I don't mean sth like that.I did write the code after a lot of effort from me(and help from others) ,but i stack in this part.I know the results because i have them in mathematica and i can see them.That's why i say the values must be around 0.3.If this post looks like you said,please the moderator deletes it. –  George Feb 22 '11 at 16:18
add comment

2 Answers

up vote 2 down vote accepted

I very strongly suggest that you try writing code without using any function handles until you're really familiar with Matlab.

The line de =@ (i,j) (2*mB*H).*mlat(i,j); is what causes your problems. In Matlab, when you define a function handle that refers to, say, an array, the function handle will use the array as it was at the time of definition. In other words, even though mlat changes inside your loop, mlat(i,j) inside the function de is always the same. In fact, you cannot even run this code unless you have previously defined mlat in the workspace.

You should therefore rewrite the main loop as follows

for iStep = 1:maxStep
   for imc = 1:mcs
      pmetro = $some function of mlat - this can be calculated using the 
                entire array as input
      %# for each element in mlat (and thus pmetro), decide whether 
      %# you have to switch the spin
      switchIdx = pmetro > 1 | pmetro < rand(size(mlat)); 
      mlat(switchIdx) = -mlat(switchIdx);
   $calculate magnetization$

Also, note that there is a command mean to take the average. No need to sum and then divide by the number of elements.

share|improve this answer
Ok,thanks Jonas!It's better to use functions,but for clarity i use function handles. –  George Feb 23 '11 at 12:37
add comment

what is mlat(i,j), Im seeing you are doing mlat(i,j)=-mlat(i,j); explain a little what you are trying to accomplish

share|improve this answer
Hello,with mlat we are creating the initial spin (L*L) lattice. (it has values -1 and 1).The mlat=-mlat is for reversing the spins.(i edited my original post,please check it) –  George Feb 22 '11 at 19:29
please describe the parameters to pmetro(parameter1,parameter2), maybe the issue is taht size, or lenghts... –  cMinor Feb 22 '11 at 19:33
At the (i,j) we store the torque at the lattice position (i,j).For the initial value we call rand and if rand<0.5 then at the atom we give torque with value 1 else -1. –  George Feb 22 '11 at 19:44
add comment

Your Answer


