最近闲来无聊,用Matlab模拟二维伊辛模型(Ising Model),Monte-Carlo,64*64个粒子,自旋有+1,-1两种取向。(后面附有代码)
照理说在无外磁场情况下,低温时所有自旋的取向应该一致,这样能量应该是最低的。
模拟得到大部分情况是这样的,可是有几个点却出现了一点点的反常现象。
下面图1是H=0时(温度,能量)曲线,可以看到低温时有几个能量凸起的情况(比如T=0.7)。然后我查看了一下T=0.7时所有粒子的自旋构型,如下图2所示,黄色+1蓝色-1,还是很有规律的。模拟得到的临界温度在T=2.3左右,与理论解析解差异不大,另外,在有外场的情况下基本不会出现这种情况,所以我对程序的正确性还是有点信心的。
然后我有以下猜想:
1)我的初始构型是随机摆放的,一半自旋向上一半自旋向下,如果初始构型全取一个方向,就不会有这个结果了(这是废话)。
2)关于一个mcs(蒙特卡洛步)的定义,一种是1mcs=遍历每个粒子做一次自旋翻转;另一种是1mcs=N次随机选取粒子做自旋翻转(N=粒子数目)。我用的是前者定义,不知各位大神推荐用哪个?
3)图2所示的状态肯定不是能量最低的状态,那它是否是亚稳态?如果是,它有何物理意义;如果不是,如何修改程序使其跑出这种构型。
附程序:
function S=MCIsing(H,T)
% Monte-Carlo simulation of 2d Ising model
% Yibing Dai 2016 April
% T=k_B*T: reduced temperature, k_B is the Boltzmann constant
% H=\mu*H: reduced magnetic field intensity, \mu is the spin magnetic moment of each atom
L=64; % length of the side of the simulation box
J=1; % exchange interaction constant
mcs=200; % Monte-Carlo sweeps
% S=zeros(L,L,mcs+1); % store the configuration of every mcs
% initial configuration: half spin up, half spin down
C=ones(L);
r=rand(L);
C(r<0.5)=-1;
S(:,:,1)=C;
% Monte-Carlo simulation
for s=1:mcs
for i=1:L
for j=1:L
DE=deltaE(C,i,j,J,H);
% determine whether or not to flip the (i,j) spin
if DE<=0
C(i,j)=-C(i,j);
else
x=rand;
if x<exp(-DE/T)
C(i,j)=-C(i,j);
end
end
end
end
S(:,:,s+1)=C;
end
function dE=deltaE(C,i,j,J,H)
% the variation of energy if flipping the (i,j) spin
L=size(C,1);
if i==1
a=C(end,j);
else
a=C(i-1,j);
end
if i==L
b=C(1,j);
else
b=C(i+1,j);
end
if j==1
c=C(i,end);
else
c=C(i,j-1);
end
if j==L
d=C(i,1);
else
d=C(i,j+1);
end
dE=2*J*C(i,j)*sum([a b c d])+2*H*C(i,j);
![二维伊辛模型有亚稳态吗?为何我模拟结果出现了亚稳态?]()
Fig1.jpg
![二维伊辛模型有亚稳态吗?为何我模拟结果出现了亚稳态?-1]()
Fig2.jpg |