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

作者: icfg66 | 来源:发表于2020-04-04 21:51 被阅读0次

Matlab蒙特卡洛模拟二维伊辛模型相变过程
科研实践的一些课题
julia - 蒙特卡洛 - 模拟布丰投针实验
Markov Chain Monte Carlo 和 Gibbs
蒙特卡洛算法随笔
蒙特卡洛理解与python实现2
Coding and Paper Letter（二十三）
二维破片撞击的Autodyn模拟
二维破片撞击的Autodyn模拟
蒙特卡洛仿真

一、什么是伊辛模型

伊辛（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格点，从同一分布开始升温，分布进行了三批温度选择：

$\beta=(0.1:0.01:0.4)$ 31个温度点， $T\in[2.5,10]$ 马尔可夫链长度为1万次，采样5万次。
$T=(0.1:0.1:2.1)$ 21个温度点，$马尔可夫链长度为1万次，采样5万次。
$T=(2.1:0.01:2.4)$ 31个温度点，马尔可夫链长度为5万次，采样25万次。

Matlab蒙特卡洛模拟二维伊辛模型相变过程
一、什么是伊辛模型伊辛（Ising）模型是描述磁系统相变最简单的模型，但模型里自旋之间的相互作用赋予了它奇妙的特...
科研实践的一些课题
一、计算二维三维伊辛模型，1、传统MATLAB；2、tensorflow能否表征临界相变？3、硬件如何优化计算？ ...
julia - 蒙特卡洛 - 模拟布丰投针实验
文章和代码参考：MATLAB模拟布丰投针实验本文适合萌新学习蒙特卡洛入门 3.0769230769230773....
Markov Chain Monte Carlo 和 Gibbs
Welcome To My Blog 一.蒙特卡洛模拟蒙特卡洛模拟(Monte Carlo Simulation...
蒙特卡洛算法随笔
蒙特卡洛蒙特卡洛可以用随机数的应用来解释，是通过随机采样来计算或模拟给定过程。经典实例就是用蒙特卡洛估计 pi ...
蒙特卡洛理解与python实现2
蒙特卡洛定义模拟是指把某一现实的或抽象的系统的某种特征或部分状态，用另一系统（称为模拟模型）来代替或模拟。为...
Coding and Paper Letter（二十三）
资源整理。 1 Coding： 1.开源项目aom，基于Fotran90编写的水光学蒙特卡洛模型。可以模拟辐射量。...
二维破片撞击的Autodyn模拟
采用Autodyn模拟一个多材料圆柱形碎片撞击板现象。主要过程打开AUTODYN 用二维轴对称建立模型求解并...
二维破片撞击的Autodyn模拟
采用Autodyn模拟一个多材料圆柱形碎片撞击板现象。主要过程打开AUTODYN 用二维轴对称建立模型求解并...
蒙特卡洛仿真
蒙特卡洛模拟(Monte Carlo Simulation)浅析 - 简书