Matlab蒙特卡洛模拟二维伊辛模型相变过程

2020-04-04 本文已影响0人 icfg66

一、什么是伊辛模型

伊辛（Ising）模型是描述磁系统相变最简单的模型，但模型里自旋之间的相互作用赋予了它奇妙的特性，最有趣的就是对称性破缺。这一模型可以被推广用于研究连续的量子相变、基本粒子超弦理论、动力学临界行为等，甚至被认为可以描述深林火灾、交通拥堵、舆论传播等社会经济现象。

伊辛模型介绍
如图，每个格点的方向只有向上或向下两者状态，但临近的自旋之间有相互作用，而且点阵可以是一维、二维、三维、甚至更高维，这两个特点让伊辛模型的严格求解成为了世纪难题。为了定量描述这个系统的能量，我们假设第个格点的自旋为，只能取+1或-1，如果相邻两个格点同方向，则它们相互作用的能量更小，设为，如果为反方向，则为，称为耦合系数，通常为正值，代表铁磁系统，如果为负值，则代表反铁磁系统。如果外磁场的强度为B，格点的自旋磁矩为，那么可以写出这个体系的哈密顿量：

$H=-J\sum_{<ij>}s_i s_j-\mu B \sum_{i=1}^N s_i$

自发对称性破缺

我们先定性地了解一下这个系统的性质，令外磁场零，当温度 $T\rightarrow 0$ 时，体系为了保持能量最低，所用的格点趋向于同方向，系统整体要么向下，要么向上，呈现强磁性。当温度 $T\rightarrow \infty$ 时，系统热运动占主导地位，格点方向呈现随机性，系统整体不带磁性，从上或从下观察体系，呈现出对称性，或者说无法通过系统磁性区分上或下。现在再考虑，当温度T从 $\infty$ 逐渐降温，那么系统必定存在某个温度 $T_c$ ，高于此温度时系统无磁性，低于此温度时，系统磁性逐渐加强。这个温度就是临界温度，也是相变点，系统从对称磁体转变为非对称磁体，而这就是对称性破缺，因为这种破缺不是外界扰动（如外加磁场）引起的，而是由内部的关联作用力造成的，所以称之为自发的对称性破缺。

自发对称性破缺
上图就是对自发对称性破缺的定性描述，当逐渐降温到时，系统磁性开始出现分化，要么向下要么向上，最终平均的磁化强度趋向于+1或-1。

我们感兴趣的问题主要有两个：第一，不同维度、不同分布的格点，其临界温度 $T_c$ 是多少；第二， $T_c$ 附近， $\overline{s}$ 随温度 $T$ 以什么样的幂指数 $\alpha$ 趋近于 $T_c$ ：
$\overline{s}\sim(T_c-T)^{\alpha}$

统计物理的思路

我们先考虑简单的9个格点的例子，实际格点数的量级为 $10^{29}$ 。

简单例子
假设格点耦合强度，那么这个9格点体系的能量为：

第一个格点有两个相邻格点，右边的与其同向，耦合能为-1，下面的与其反向，耦合能为1；第二个格点有三个相邻格点，左和下与其同向，耦合能为-2，右边与其反向，耦合能为1。类似的可以计算其余格点的耦合能。最后除以2是因为每个相互作用重复计算了一次。而这一特定分布以某概率P出现：

对于9个格点的体系，总共有种分布，每种分布的出现概率于其总能量有关，所以分布概率满足归一化条件。给定不同温度，我们可以计算出不同温度下平均磁化率的数学期望，得到的曲线。

但是对于粒子为 $10^{29}$ 量级的格点，可能的分布有 $2^{10^{29}}$ 种，根本不可能统计出结果。

一维的情况可以通过数学上的处理，最终可以提出N，并得到 $T_c=0$ ，也就是说一维的伊辛模型不会有自发的对称性破缺，这是因为一维的格点只有两个相邻格点，相互作用太弱，不足以对抗热运动，始终表现为整体0磁化率的对称状态。

二维的情况下，如果用平均场近似的方法（具体可以参考林宗涵的热力学与统计物理），基本思想是将相互作用的耦合能转化为外磁场强度，这就可以用近独立的模型来计算配分函数，进而得到所有的统计量，获得的临界温度为 $T_c= \frac{2J}{k}$ ，平均磁化率：
$\overline{s}\sim(T_c-T)^{1/2} ~~~~~~~~(T\rightarrow T_c^-)$

1944年，昂萨格推导出了二维伊辛模型的严格解，临界温度 $T_c=\frac {2.269J}{k}$ ，平均磁化率：
$\overline{s}\sim(T_c-T)^{1/8} ~~~~~~~~~(T\rightarrow T_c^-)$

二、二维伊辛模型的精确解

平均磁化率 $\overline s$

这里只列出二维伊辛模型精确解的结论，推导过于复杂，定义 $\beta\equiv \frac{1}{kT}$ ，则平均磁化率：
$\overline{s} =\begin{cases} 0,&{T>T_c} \\ [1-\sinh(2\beta J)^{-4}]^{1/8}, &{T\leq T_c}\end{cases}$
令 $\sinh(2\beta J)=0$ 可以解得： $T_c=\frac{2J}{k\ln(1+\sqrt{2})}\approx \frac{2.269J}{k}$ ，令 $T=T_c-\delta T$ ，小量泰勒展开化简可以得到：
$\overline{s}=[4·\frac{2J}{kT_c}\cosh(\frac{2J}{kT_c})·\frac{\delta T}{T_c}]^{1/8}=1.224[1-\frac{T}{T_c}]^{1/8}\sim(T_c-T)^{1/8}$
当 $T_c \rightarrow 0$ 时，容易得到结果为1，表面所有的格点同向，平均磁化率为1。

假设耦合强度 $J$ 等于1开尔文温度下的热运动能量，即 $J=k$ ，做出 $\overline{s}-T$ 曲线如下：

平均磁化强度s-T

平均能量和比热

另一个能够判断相变的参数是比热 $C_v=\frac{d \overline{E}}{dT}$ ，这里的 $\overline{E}$ 表示单个格点的平均能量，定性来看，当温度趋于0时，所有格点同向， $\overline{E}=-4/2=-2$ ，当温度趋于无穷时，格点方向随机，某格点四周平均有两个同向和两个方向， $\overline{E}=0$ 。其曲线如下(令 $k=J$ )：

能量与比热随温度的变化图

这里的比热在临界温度时会突然增大，表面临界温度附近变化有个小温度变化，需要吸收极大能量，这也很符合相变的特点。具体的计算公式和matlab代码详见附录A。

虽然大体了解了相变的过程，以及理论上的精确解，我们能否通过实验的方式，验证这一结论呢？借助计算机模拟这一过程来验证结果呢？

三、二维伊辛模型模拟

因为不可能遍历所有的格点组合，我们只能利用采样的方式去计算平均能量，采样的条件应该是体系在某个温度下已经平衡。计算机模拟的基本思想是，首先随机给定一种分布，在特定温度下，让体系趋向平衡，再在这个平衡体系中采样求平均。

假设体系有 $20\times 20$ 的格点，初始时同一分布，相当于温度很低；
我们设定一个希望“加热”到的平衡温度 $T_0$ ，接下来是模拟最关键的地方，如何改变格点的分布以趋向于设定温度 $T_0$ ？
1、我们先任意选择一个格点，计算改变这个格点的能量变化 $\Delta E$ ，因为体系出现的概率正比于 $e^{-E/kT}$ ，那么格点变化前后两个体系出现的概率比为 $e^{-\Delta E/kT}$ ，或者说格点改变的概率与不改变的概率比为： $1:e^{-\Delta E/kT}$ ，那么格点需要改变的概率为 $e^{-\Delta E/kT}$ ，因此我们产生一个(0,1)概率随机数，如果它小于 $e^{-\Delta E/kT}$ 则选择改变格点；
2、重复1的步骤，直到体系相对平稳，这个过程称为马尔可夫链，之后在平稳的体系下采样若干次并做统计平均，获得平均能量 $\overline E$ ，平均磁化率 $\overline s$ ，以及描述能量方差的比热 $C_v$ 。
3、设定另一个希望考察的温度，重复1、2步骤，获得该温度下的统计参数，并和理论值比较。

我们同样假设 $J=k$ ，选取格点数为 $20\times 20$ 。临界温度点附近，马尔可夫链长取5万次，采样数为25万次；其他温度点马尔可夫链长1万次，采样数为5万次。这是因为临界温度附近的涨落很大，需要更长的时间趋向平衡，需要更多的统计样本获得较准确平均值。详细的代码及解释可以参看附录B。

模拟平均磁化率 $\overline{s}-T$

平均磁化率的模拟结果

可以看出平均磁化率在临界温度附近很不稳定，这是因为临界相变时涨落很大的缘故，高温时的磁化率不是严格的零，可能与格点数少和马尔可夫链较短有关系，如何确定 $T_c$ 呢？通过平均磁化率求 $T_c$ 比较困难，一般是通过比热 $C_v$ 发散的位置确定 $T_c\approx 2.3$ ，参考下一部分。

选取T=1.7~2.2的16个数据点，拟合曲线：
$\ln\overline{s}\sim \alpha \ln(T_c-T)$
得到 $\alpha \approx 0.12,R^2=0.89$ ，这和理论值 $1/8=0.125$ 相当接近。

模拟平均能量 $\overline{E}-T$

平均能量和比热的模拟结果

在临界温度附近进行了较密集的温度取点，而且加长了马尔可夫链，但是仍能看到较大的涨落，可以通过这个现象来确定临界温度 $T_c\approx 2.3$ 。

最后，让我们欣赏一下格点从同一分布到临界温度的变化过程吧，为了便于观察，选择40X40的格点，马尔可夫链长为1万：

初始同一分布的格点，逐渐趋向于高温 $T=5$ 下的平衡态，格点最后呈现随机分布：

T=5，初始为有序

初始同一分布的格点，温度降低至 $T=3$ 的平衡态，格点最后呈现小块状：

T=3，初始为有序

初始无序分布的格点，温度降低到临界温度 $T=2.3$ 时，格点最后呈现的块状增大。

T=2.3，初始为无序

初始无序分布的格点，温度降低到临界温度以下 $T=2$ ，格点最后呈现的大的块状，说明已经发生了明显的相变。

T=2，初始为无序

初始无序分布的格点，温度降低更多至 $T=1$ ，格点越来越趋向于同一分布。

T=1，初始为无序

附录

A、平均能量和比热的精确解：(前方高能)

定义 $\beta \equiv 1/kT$
$\overline E=-J\cot(2\beta J)\times[1+\frac{2}{\pi}A·B(\lambda)]$
$\begin{cases} A\equiv 2\tanh(2\beta J)^2-1\\ B(\lambda)\equiv \int_0^{\pi/2} \frac{d\phi}{\sqrt{1-\lambda^2\sin\phi ^2}}\\ \lambda\equiv \frac{2\sinh(2\beta J)}{\cosh(2\beta J)^2} \end{cases}$
帝国主义都是纸老虎，我们仔细发现，只要确定了温度 $T$ 或 $\beta$ ，那么可以依次确定 $\lambda,B(\lambda),A,\overline E$ ，也就是平均能量和温度是一一对应的，最后通过求导得到比热 $C_v=d \overline E/d T$ 。

%matlab code
clear;clc;
beta=0.1:0.01:1; %温度从10到1
for i=1:1:size(beta,2);%遍历beta
lambda=2.*sinh(2.*beta(i))./cosh(2.*beta(i)).^2;%计算lambda
phi=linspace(0,pi/2,1000);%求B的积分参数
b=1./sqrt(1-lambda.^2.*sin(phi).^2);
B=trapz(phi,b);%积分B(lambda)
A=2*tanh(2.*beta(i)).^2-1;
e_bar(i)=-coth(2.*beta(i)).*(1+2/pi.*A.*B); %每个格点的平均能量
end
%%
plot(1./beta,e_bar,'k','LineWidth',2);hold on;%E-T曲线
beta1=beta(2:end)/2+beta(1:end-1)/2;
cv=-beta1.^2.*diff(e_bar)./diff(beta); %对能量求导得到比热
plot(1./beta1,cv,'r','LineWidth',2);

平均能量和比热随温度的变化图

B、蒙特卡洛马尔可夫链模拟二维伊辛模型相变过程

%matlab code
clear;clc;
n=10000;                       %马尔可夫链长度1万
ns=20;                          %20*20的格点 
beta_mc=(0.1:0.01:0.4);         %温度从10到2.5,链长1万，样品长5万
%T_mc=(2.1:0.01:2.4);          %第二次模拟温度设定，临界温度附近取点更密集，还要调整n=50000
%beta_mc=1./T_mc;
tic;                               %计时用，n=10000时，通常需要跑一分多钟
for jj=1:1:size(beta_mc,2)
X=sign(rand(ns,ns));        %方向一致，相当于从0度开始升温
%马尔可夫链长度为5万次
for j=1:1:n
    %随机选取一个格点，行列存储在index[1,2]
    index=unidrnd(ns,1,2);      
    % 利用周期性边界条件，分别计算格点上下左右四个点行列坐标
    tmp1=rem(index(1),ns)+1;tmp2=rem(index(1)+1,ns)+1;tmp3=rem(index(1)-1,ns)+1;
    tmp4=rem(index(2),ns)+1;tmp5=rem(index(2)+1,ns)+1;tmp6=rem(index(2)-1,ns)+1;
    % 计算改变格点方向后的能量变化
    cen=X(tmp1,tmp4);right=X(tmp1,tmp5);left=X(tmp1,tmp6);
    up= X(tmp2,tmp4);down= X(tmp3,tmp4);
    deE=2*cen*(right+left+up+down);
    % 判断是否改变格点
    if rand<exp(-deE*beta_mc(jj))
        X(tmp1,tmp4)=-X(tmp1,tmp4);
    end
end    

% 取样5万次，平衡时同样需要判断是否改变格点
for j=1:1:5*n
    index=unidrnd(ns,1,2);
    % 利用周期性边界条件，分别计算格点上下左右四个点行列坐标
    tmp1=rem(index(1),ns)+1;tmp2=rem(index(1)+1,ns)+1;tmp3=rem(index(1)-1,ns)+1;
    tmp4=rem(index(2),ns)+1;tmp5=rem(index(2)+1,ns)+1;tmp6=rem(index(2)-1,ns)+1;
    % 计算改变格点方向后的能量变化
    cen=X(tmp1,tmp4);right=X(tmp1,tmp5);left=X(tmp1,tmp6);
    up= X(tmp2,tmp4);down= X(tmp3,tmp4);
    deE=2*cen*(right+left+up+down);
    % 判断是否改变格点
    if rand<exp(-deE*beta_mc(jj))
        X(tmp1,tmp4)=-X(tmp1,tmp4);
    end
    %计算一种特定分布时的平均磁化率
    m(j)=abs(mean(mean(X)));  
    %计算一种特定分布时的平均能量
    Xt1=X;Xt1(1,:)=[];Xt1=[Xt1; X(1,:)];
    Xt2=X;Xt2(:,1)=[];Xt2=[Xt2, X(:,1)];
    e(j)=-mean(mean(X.*Xt1+X.*Xt2));
    e2(j)=e(j)^2;
end
% 特定温度下的统计量
m_bar(jj)=mean(m);   
e_bar(jj)=mean(e);
cv_bar(jj)=beta_mc(jj)^2*ns^2*std(e)^2;
end
toc;
% 作图观察
figure(1);
plot(1./beta_mc,m_bar,'ko');
figure(2);
plot(1./beta_mc,e_bar,'ko');
figure(3);
plot(1./beta_mc,cv_bar,'ro');

20X20格点，从同一分布开始升温，分布进行了三批温度选择：