ising模型c语言程序,Ising 模型

本文深入探讨了Ising模型,一种用于描述物质铁磁性的数学模型。涵盖了从一维到二维模型的计算方法,包括Weiss平均场理论和Monte Carlo模拟,并提供了具体的计算实例。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

66b52468c121889b900d4956032f1009.png

8种机械键盘轴体对比

本人程序员,要买一个写代码的键盘,请问红轴和茶轴怎么选?

Ising模型是一个以物理学家恩斯特·易辛为名的数学模型,用于描述物质的铁磁性。一个二维的方晶格易辛模型是已知最简单而会产生相变的物理系统之一。

模型描述:

N个经典粒子(编号),其哈密顿量只与每个粒子的自旋自由度(只有单方向分量)有关,与空间运动自由度()无关;粒子静止的摆放在周期性d维晶格上(冻结了空间运动自由度)。

体系微观状态的描述为: 可见体系共有种微观状态(相空间的大小)

Ising Hamiltonian:

<>代表不重复的最近邻的键 ,最近邻格点上的自旋与自旋之间相互作用,自旋与外场之间有耦合。J>0称为铁磁Ising模型,J<0称为反铁磁Ising模型。

一维Ising模型计算

一维Ising模型定性理解

零场零温(T=0,B=0)下

T=0,基态二重简并

aGc11ciAmI.png

T>0,基态与激发态都以一定概率出现:

cJLLIjLGGL.png

温度会打破完全有序的状态,给体系带来熵。第一激发态N-1重简并。

一般来说,自由能F=U-TS取极小决定了体系的热力学状态,这里面有序项,无序项(只考虑元激发)。相变点即,可以得到,

热力学极限下,,即一维Ising系统不存在非零温相变。

一维Ising模型严格解

具体看体系热力学性质需要严格求出体系配分函数,这里需要用到转移矩阵的方法,哈利顿量写成

配分函数

于是可以写为:

其中称为转移矩阵,为温度相关系数。又已知,任意矩阵 所以得到

其中 通过其得到体系热力学性质。

热力学极限下,可以得到

自由能是温度T与磁场B的解析函数,无奇异性,任意阶可微,不存在自发磁化。

二维Ising模型计算

Weiss平均场理论计算

除了少数几个特殊情况以外,无法得到Ising模型的严格解。忽略格点自旋的相关性。 通过配分函数计算平均磁矩,得到自洽条件, 考虑自发磁化()

由自洽条件得到居里温度

Ising模型经典计算方法

Metropolis算法

Metropolis–Hastings 算法是在数值模拟易辛模型时最常用的一种蒙特卡洛方法。首先取一个冻结自旋系统模型,初始构型的选择依赖于温度。高温下可以选择完全随机无序的初始构型,低温下选择有序铁磁或者反铁磁构型。

随机或逐个选取一个或几个自旋进行反转形成新构型,将一个老构型改变成另外一个新构型时,新构型被认可概率正比于,其中为两构型能量差。

如果新构型能量低,接受该构型。如果新构型能量较高,则产生一个0~1随机数,如就接受,否则维持原构型。

核心子程序代码

function m=Isingrand(T)

N=1000;

s=zeros(N,N);

%T is temperature

J=1;%exchange interaction strength

k=1;%Boltzmann constant

for i=1:N

for j=1:N

if(rand>0.5)

s(i,j)=1;

else

s(i,j)=-1;

end

end

end

for i=1:2000*N^2

[a,x,y]=spinflip(s,k,T,N,J);

s(x,y)=s(x,y)*a;

end

m=abs(sum(sum(s))/N^2);

end

function [a,x,y]=spinflip(s,k,T,N,J)

x=floor(rand*N+1);

y=floor(rand*N+1);

if x~=1&&x~=N&&y~=1&&y~=N

dH=2*J*s(x,y)*(s(x-1,y)+s(x+1,y)+s(x,y-1)+s(x,y+1));

elseif x==1&&y~=1&&y~=N

dH=2*J*s(1,y)*(s(1,y-1)+s(1,1+y)+s(2,y)+s(N,y));

elseif x==N&&y~=1&&y~=N

dH=2*J*s(N,y)*(s(N,y-1)+s(N,1+y)+s(1,y)+s(N-1,y));

elseif x~=1&&x~=N&&y==1

dH=2*J*s(x,1)*(s(x-1,1)+s(x+1,1)+s(x,2)+s(x,N));

elseif x~=1&&x~=N&&y==N

dH=2*J*s(x,N)*(s(x-1,N)+s(x+1,N)+s(x,1)+s(x,N-1));

elseif x==1&&y==1

dH=2*J*s(1,1)*(s(2,1)+s(N,1)+s(1,2)+s(1,N));

elseif x==1&&y==N

dH=2*J*s(1,N)*(s(2,N)+s(N,N)+s(1,1)+s(1,N-1));

elseif x==N&&y==N

dH=2*J*s(N,N)*(s(1,N)+s(N-1,N)+s(N,1)+s(N,N-1));

elseif x==N&&y==1

dH=2*J*s(N,1)*(s(1,1)+s(N-1,1)+s(N,2)+s(N,N));

end

a=1;

if(dH<=0)

a=-1;

elseif (exp(-1*dH/k/T)>rand)

a=-1;

end

end

计算结果

由于计算模型较大,自己计算机运行能力有限,没有执行太多循环。通过两张图,定性的可以开出与的正相关关系。

左图,右图,均为磁化率-温度()图像。

jF3D0dCalK.png

Ising模型量子模拟方法

二维伊辛模型是一个重要的统计物理模型,用于描述强关联体系的基态和激发态。在无外加磁场下,其哈密顿量可以表示为: $$ H = -J\sum_{<i,j>}s_is_j $$ 其中 $s_i$ 表示在位置 $i$ 的自旋,取值为 $\pm 1$;$<i,j>$ 表示 $i$ 和 $j$ 为相邻格点,$J$ 是自旋之间的交换作用强度。 蒙特卡罗方法是一种常用的计算机算法,可以模拟随机过程,因此非常适合用于二维伊辛模型的模拟。下面我们基于 Matlab 实现无外加磁场下的二维伊辛模型的模拟,并绘制能量、磁化强度、热容和磁化率随度的变化。 首先,我们需要初始化模型。我们选取一个 $L\times L$ 的正方形格点,每个格点上有一个自旋,初始状态可以随机生成。为了方便计算,我们可以使用周期性边界条件,即相邻的格点相互连接,形成一个环。 ```matlab % 初始化二维伊辛模型 L = 32; % 格点数 T = linspace(1.5,3,30); % 度范围 J = 1; % 交换作用强度 N = L*L; % 自旋数 s = sign(rand(N,1)-0.5); % 随机生成初始状态 ``` 接下来,我们需要使用 Metropolis 算法进行模拟。具体地,我们每次随机选取一个自旋,计算翻转该自旋后的能量变化 $\Delta E$。如果 $\Delta E<0$,则接受该翻转;否则以一定概率 $\exp(-\Delta E/T)$ 接受该翻转。这里的 $T$ 是度,控制了模型的热力学行为。 ```matlab % Metropolis 算法模拟二维伊辛模型 ntherm = 1000; % 热化步数 nmc = 10000; % 蒙特卡罗步数 E = zeros(length(T),nmc); % 能量 M = zeros(length(T),nmc); % 磁化强度 C = zeros(length(T),1); % 热容 X = zeros(length(T),1); % 磁化率 for k=1:length(T) for i=1:ntherm r = randi([1,N]); % 随机选择一个自旋 s0 = s; s(r) = -s(r); % 翻转该自旋 deltaE = 2*J*s(r)*(s(mod(r-2,N)+1)+s(mod(r,N)+1)+s(mod(r,L)+r-L)+s(mod(r-1,L)+r-1-L)); if deltaE > 0 && rand > exp(-deltaE/T(k)) s(r) = s0(r); % 恢复原状态 end end E(k,1) = -J*sum(s.*(s(mod(1:N-1,N)+1)+s(mod(1:N-L-1,N)+L+1))); M(k,1) = sum(s)/N; for i=2:nmc r = randi([1,N]); % 随机选择一个自旋 s0 = s; s(r) = -s(r); % 翻转该自旋 deltaE = 2*J*s(r)*(s(mod(r-2,N)+1)+s(mod(r,N)+1)+s(mod(r,L)+r-L)+s(mod(r-1,L)+r-1-L)); if deltaE > 0 && rand > exp(-deltaE/T(k)) s(r) = s0(r); % 恢复原状态 end E(k,i) = E(k,i-1) + deltaE; M(k,i) = M(k,i-1) + 2*s(r)/N; end E2 = E(k,:).^2; M2 = M(k,:).^2; C(k) = (mean(E2)-mean(E(k,:))^2)/(T(k)^2); X(k) = (mean(M2)-mean(M(k,:))^2)/T(k); end ``` 最后,我们可以绘制能量、磁化强度、热容和磁化率随度的变化曲线。这里我们使用 Matlba 的 `plot` 函数进行绘图。 ```matlab % 绘制能量、磁化强度、热容和磁化率随度的变化曲线 figure; subplot(2,2,1); plot(T,E/N,'.-'); xlabel('Temperature'); ylabel('Energy per spin'); subplot(2,2,2); plot(T,abs(M),'.-'); xlabel('Temperature'); ylabel('Magnetization per spin'); subplot(2,2,3); plot(T,C,'.-'); xlabel('Temperature'); ylabel('Specific heat'); subplot(2,2,4); plot(T,X,'.-'); xlabel('Temperature'); ylabel('Susceptibility'); ``` 运行以上代码,即可得到能量、磁化强度、热容和磁化率随度的变化曲线。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值