ising model from mathematica to matlab

Question

i have this problem from mathematica,which refers to ising model.Here is a piece of code:

J1k = Table[2 RandomInteger[] - 1, {L}, {L}];
J2k = Table[2 RandomInteger[] - 1, {L}, {L}];

energy1 := 
 Module[{ii1, ii2, jj1, jj2, kk1, kk2, dG1, ener1, ener2, ener}, 
  ener = 0.;
  Do[
   Do[
    jj1 = ii1 + 1; jj2 = ii2 + 1; kk1 = ii1 - 1; kk2 = ii2 - 1;
    If[jj1 > L, jj1 = 1]; If[jj2 > L, jj2 = 1];
    If[kk1 < 1, kk1 = L]; If[kk2 < 1, kk2 = L];
    ener1 = -J1k[[ii1, ii2]]*mlat[[ii1, ii2]]*mlat[[jj1, ii2]] - 
      J1k[[kk1, ii2]]*mlat[[ii1, ii2]]*mlat[[kk1, ii2]];
    ener2 = -J2k[[ii1, ii2]]*mlat[[ii1, ii2]]*mlat[[ii1, jj2]] - 
      J2k[[ii1, kk2]]*mlat[[ii1, ii2]]*mlat[[ii1, kk2]];
    ener = ener1 + ener2 - 2*mlat[[ii1, ii2]] *\[Mu]Bk,
    {ii1, 1, L}],
   {ii2, 1, L}];
  ener = ener/2.]

energy = energy1

Here is what i have done:

J1k=2*randint(L,L)-1;
J2k=2*randint(L,L)-1;
I have created a function:
function Energy1 =energy1()
%spin interaction with the neighbors
ener=0;
L=16;
J1k=2*randint(L,L)-1;
J2k=2*randint(L,L)-1;
for ii2=1:L
for ii1=1:L
    jj1=ii1+1;jj2=ii2+1;kk1=ii1-1;kk2=ii2-1;
    if (jj1>L & jj1==1)
    end
        if (jj2>L & jj2==1)
        end
            if (kk1<1 & kk1==L)
            end
                if (kk2<1 & kk2==L)
                end
                    ener1=-J1k(ii1,ii2).*mlat(ii1,ii2).*mlat(jj1,ii2)-J1k(kk1,ii2).*mlat(ii1,ii2).*mlat(kk1,ii2);
                    ener2=-J2k(ii1,ii2).*mlat(ii1,ii2).*mlat(ii1,jj2)-J2k(ii1,kk2).*mlat(ii1,ii2).*mlat(ii1,kk2);
                    ener=ener1+ener2-2.*mlat(ii1,ii2).*mBk;


end                
end
ener=ener ./2;

end

energy=@energy1

The problem is that when i call energy "energy=@energy1" no value returns.It supposes to return a value "1" or "-1".

Also, it continues (in mathematica ) :

T = 5.2; dT = 0.1; Umean = {};
energy = energy1;
Do[energyseries = {}; T = T - dT;
  Do[i1 = RandomInteger[{1, L}]; i2 = RandomInteger[{1, L}];
   j1 = i1 + 1; j2 = i2 + 1; k1 = i1 - 1; k2 = i2 - 1;
   If[j1 > L, j1 = 1]; If[j2 > L, j2 = 1];
   If[k1 < 1, k1 = L]; If[k2 < 1, k2 = L];
   dG1 = -J1k[[i1, i2]]*mlat[[j1, i2]] - 
     J1k[[k1, i2]]*mlat[[k1, i2]];
   dG2 = -J2k[[i1, i2]]*mlat[[i1, j2]] - 
     J2k[[i1, k2]]*mlat[[i1, k2]];
   dE = 2.*mlat[[i1, i2]] (dG1 + dG2 + \[Mu]Bk);
   W = N[Exp[-dE/T]];
   If[W < 1 && W > Random[] || dE <= 0, 
    mlat[[i1, i2]] = -mlat[[i1, i2]]; index = 1, index = 0];
   energy = energy + dE*index // N;
   AppendTo[energyseries, energy],
   {i, 1, Maxstep}];
  PrependTo[Umean, {T, Mean[energyseries]/L2}],
  {jT, 1, 51}];

And i did :

T=5.2;
dT=0.1;
Umean=zeros();
energy=@energy1
for jT=1:51
    energyseries=zeros();
    T=T-dT;
    for i=1:Maxstep
    i1=randint(1,L);
    i2=randint(1,L);
    j1=i1+1;
    j2=i2+1;
    k1=i1-1;
    k2=i2-1;
    if (j1>L & j1==1)
    end
        if (j2>L & j2==1)
        end
            if (k1<1 & k1==L)
            end
                if (k2<1 & k2==L)
                end
                    dG1=-J1k(i1,i2).*mlat(j1,i2)-J1k(k1,i2).*mlat(k1,i2);
                    dG2=-J2k(i1,i2).*mlat(i1,j2)-J2k(i1,k2).*mlat(i1,k2);
                    dE=2.*mlat(i1,i2).*(dG1+dG2+mBk);
                    W=exp(-dE./T);
                    if (W<1 & W>rand() | dE<=0)
                        mlat(i1,i2)=-mlat(i1,i2);
                        index=1;
                        index=0;
                    end
                    energy=energy+dE.*index;
                    energyseries(:)=[energy]
    end
    Umean(:)=[T mean(energyseries(:))./L2]
end

which gives me message " Subscript indices must either be real positive integers or logicals.

Error in ==> ising at 58 dG1=-J1k(i1,i2).*mlat(j1,i2)-J1k(k1,i2).*mlat(k1,i2);" Any ideas?

Ghaul · Accepted Answer · 2011-02-19T17:01:57.490

1

For your first problem, you haven't set a value for the output, "Energy1". Just do

Energy1 = ener;

at the end of the function.

I think that your second problem is your use of the "if" statements. "If" in Matlab does not work like "if" in mathematica. In Matlab it works like this:

if (expression==true)   
     execute statement
end

So what you want to do is

if (j1>L)
     j1 = 1;
end

etc...

EDIT to summarize discussion in comments: Use "rand(L)" to define i1 and i2 instead of "randint"

edited Feb 19 '11 at 17:01

answered Feb 19 '11 at 14:07

Ghaul

3,340
1
19
24

For the first problem ,the "Energy1=ener"doesn't work too.It gives me "energy=@energy1".For the second problem you are right!It seems i was totally confused!I edited the loops as you said,but it still gives me "Subscript indices must either be real positive integers or logicals...Error in ==> ising at 62 dG1=-J1k(i1,i2).*mlat(j1,i2)-J1k(k1,i2).*mlat(k1,i2); – George Feb 19 '11 at 14:25
I think you use the "@" wrong. When you write "energy = @energy1" you define "energy" to be the function "energy1". You can then call upon it by writing "energy()". To call upon the function "energy1" directly, just write "energy = energy1()". The index error is because Matlab starts indexing at 1 and not 0. Your values k1 and k2 can take the value 0. You then get an error when you try to access J1k in index 0. – Ghaul Feb 19 '11 at 14:57
If i do "energy=@energy1"and then "energy()" it gives me the same as before "enery=@energy1".If i do "energy=energy1()" it doesn't recognize any of the variables i have in the function file(as L,mlat).As for the index error which you said depends from the values k1,k2 how can i overcome that? – George Feb 19 '11 at 15:12
Functions in Matlab can't access the variables in your workspace. You either need to define them inside the function or give them as input. I suggest you edit the line where your define your function and write "function Energy1 = energy1(L,mlat)" and then call on it with "energy=energy1(L,mlat);" For the indexing problem, I'll have to think about it a bit.. – Ghaul Feb 19 '11 at 15:29
The problem with energy1 ok ,now it works.As for the values(all have L=16 values) :i1 and i2 takes 0 and 1 , j1 and j2 takes 1 and 2 and k1 and k2 takes 0 and -1 – George Feb 19 '11 at 15:49
Good. The problem is in the function randint, which doesn't work like Mathematicas RandomInteger. You should use "randi(L)" instead of "randint" – Ghaul Feb 19 '11 at 15:50
I have used ranint until now in some problems and seems to be ok..Anyway,i used randi but now the i1,i2,j1,j2,k1,k2 have LxL lines and columns and also it gives me the message "??? Error using ==> minus Matrix dimensions must agree. Error in ==> ising at 63 dG1=-J1k(i1,i2).*mlat(j1,i2)-J1k(k1,i2).*mlat(k1,i2); " – George Feb 19 '11 at 16:02
L=16 right? Then randint(1,16) will give you a vector of length 16 with only 0s and 1s. Writing randi(16) gives you a random number between 1 and 16, which I think is what you want? It could be that you have an older version of Matlab where randint works differently from mine.. – Ghaul Feb 19 '11 at 16:11
If i have randint(1,L) the equivalent of randi is randi(L)?And if i have randint(L,L) the randi will be?As for the error it says"minus matrix dimensions must agree" is from not defining randi properly? – George Feb 19 '11 at 16:17
I don't understand. According to the Mathematica code, you want i1 and i2 to be a random number between 1 and 16? You do this be setting "i1=randi(16)" (not randi(1,16),which will give a 16x16 matrix, which I suspect is the reason for your error.) – Ghaul Feb 19 '11 at 16:28
For defining J1k and J2k it is correct to use randint, since RandomInteger[] gives 0 or 1, but for i1 and i2 you have to use randi since RandomInteger[{1,16}] gives a random number between 1 and 16. – Ghaul Feb 19 '11 at 16:37
Ok ,i got it.I will try to fix it now!Thanks for the help!Also,can you tell me how can i do this from mathematica in matlab? "T = Table[{Null, Null}, {Maxstep}]"If i do "T=NaN(1,Maxstep)"it will give me 1 line and nan values ,but i want pairs of nan values. – George Feb 19 '11 at 16:45
Also,whenever you can please check http://stackoverflow.com/questions/5032029/map-command-and-few-other-from-mathematica-to-matlab ,i will be grateful...I swear it will not take you time! :) – George Feb 19 '11 at 16:46
Great=) Do you mean like this: T = NaN(1,2,Maxstep) or this: T = NaN(2,Maxstep)? – Ghaul Feb 19 '11 at 16:59
I don't have time now, but might check it out tomorrow. – Ghaul Feb 19 '11 at 17:00
I think NaN(1,2,Maxstep) is what i want.I'll see how it goes.Thank you very much for all your help! – George Feb 19 '11 at 17:15

nibot · Answer 2 · 2011-02-19T22:06:18.820

I implemented an Ising model in Matlab a few years ago. Perhaps my code would be useful to you; it is available in this note: Monte Carlo investigation of the Ising model (PDF). The code is very short and Matlabesque and begins on page 6.

As an example of idiomatic Matlab programming, suppose you have a matrix grid which encodes your grid of spins as either +1 (up) or -1 (down). Then you can compute the average magnetization of every site's neighbors without using any for-loops by doing:

% Calculate the number of neighbors of each cell
neighbors = circshift(grid, [ 0 1]) + ...
            circshift(grid, [ 0 -1]) + ...
            circshift(grid, [ 1 0]) + ...
            circshift(grid, [-1 0]);

Then you can compute the change in energy due to flipping any spins:

% Calculate the change in energy of flipping a spin
DeltaE = 2 * J * (grid .* neighbors);

This statement works on every grid site simultaneously; such "vectorization" (computing over entire arrays at once) is the key to efficient Matlab code.

ising model from mathematica to matlab

2 Answers2