Zweidimensionales Ising-Modell

© Prof. Dr. Jürgen Schnack, Herbstakademie 2007

Ising-Modell auf einem Quadratgitter der Größe NMax * NMax.
Die Temperatur kann vorgegeben werden.
Die Wechselwirkung der Spins ist 1.

Einstellungen

NMax=10;
T=2.6;

Wechselwirkungsmatrix

NMaxSq=(NMax^2);

Spins sind wie folgt abgezählt:
(({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}))

Hilfsfeld, in dem die Spinnummer in den Platz umgerechnet wird

Platz=Table[{IntegerPart[((k-1)/NMax)]+1,Mod[k-1,NMax]+1},{k,1,NMaxSq}]

Wechselwirkungsmatrix JMat

JMat=Table[0,{i,1,NMaxSq},{j,1,NMaxSq}];
For[i1=1,i1<NMax+1,i1++,
For[j1=1,j1<NMax+1,j1++,
For[i2=1,i2<NMax+1,i2++,
For[j2=1,j2<NMax+1,j2++,
(* Print[i1," ",i2," ",j1," ",j2," ",(i1-1)*NMax+j1," ",(i2-1)*NMax+j2]; *)
If[Mod[Abs[i1-i2],NMax]==1&&j1==j2,JMat[[(i1-1)*NMax+j1,(i2-1)*NMax+j2]]=1];
If[Mod[Abs[i1-i2],NMax]==NMax-1&&j1==j2,JMat[[(i1-1)*NMax+j1,(i2-1)*NMax+j2]]=1];
If[Mod[Abs[j1-j2],NMax]==1&&i1==i2,JMat[[(i1-1)*NMax+j1,(i2-1)*NMax+j2]]=1];
If[Mod[Abs[j1-j2],NMax]==NMax-1&&i1==i2,JMat[[(i1-1)*NMax+j1,(i2-1)*NMax+j2]]=1]
]
]
]
];
MatrixForm[JMat];

Ausgangszustand

Spins=Table[2 RandomInteger[]-1,{i,1,NMax},{j,1,NMax}];
Spins0=Spins;
Spins

I1=Show[Graphics[Raster[Spins]/.{1->{1,0,0},-1->{0,0,1}},AspectRatio->1]]

Export["I1.jpg",I1];

Metropolis

MaxSweeps=100;

DeltaE=(ENeu-EAlt) ist zweimal die Summe, weil man -s_i*Summe-s_i*Summe hat.

Spins=Spins0;
For[Sweeps=1,Sweeps<MaxSweeps+1,Sweeps++,For[i=1,i<NMaxSq+1,i++,DeltaE=2 Spins[[Platz[[i,1]],Platz[[i,2]]]] (∑_(k=1)^NMaxSq)JMat[[i,k]] Spins[[Platz[[k,1]],Platz[[k,2]]]];If[DeltaE≤0,Spins[[Platz[[i,1]],Platz[[i,2]]]]=-Spins[[Platz[[i,1]],Platz[[i,2]]]],If[RandomReal[]≤(E^(-(DeltaE/T))),Spins[[Platz[[i,1]],Platz[[i,2]]]]=-Spins[[Platz[[i,1]],Platz[[i,2]]]]]]];
Print[Graphics[Raster[Spins]/.{1->{1,0,0},-1->{0,0,1}},AspectRatio->1]]];