1、热热 力力 学与学与 统计统计 物理物理 课课 程程 设计设计计算机模拟计算机模拟(Matlab 实现实现 )题目: Ising 模型铁磁-顺磁相变现象的 Monte Carlo 模拟学院: 数理学院 专业: 06 应用物理 指导老师: 刘旭东 同组人员: 姓 名 学 号彭加福 0640502112黎 婕 0640502102胡进军 0640502107热力学统计物理计算机模拟- 2 -Ising 模型铁磁顺磁相变现象的 Monte Carlo 模拟 （江苏科技大大学 06 应用物理）摘要：课程设计中首先采用 Ising model 的思想建立一个二维的模型，然后利用重要性抽样和 Monte
2、Carlo方法及其思想模拟铁磁-顺磁相变过程。计算了顺磁物质的能量平均值 Ev、热容 Cv、磁化强度 M 及磁化率X 的值，进而研究 Ev、Cv、M、X 与温度 T 的变化关系并绘制成 Ev-T 图、Cv-T 图、M-T 图、X-T 图，得出顺磁物质的内能随着温度的升高先增大而后趋于稳定值；热容 Cv、磁化率 X 随着温度的升高先增大后减小；磁化强度 M 在转变温度 Tc 处迅速减小为零，找出铁磁相变的转变温度 Tc 大约为 2.35。关键字：蒙特卡罗方法、伊辛模型、重要性抽样、顺磁相变、居里点、磁化率、磁 畴一、内容概述1907 年 法 国 科 学 家 外 斯 系 统 地 提 出 了 铁 磁
3、 性 假 说 ， 其 主 要 内 容 有 ： 铁 磁 物 质 内 部 存 在 很 强 的“分 子 场 ”， 在 “分 子 场 ”的 作 用 下 ， 原 子 磁 矩 趋 于 同 向 平 行 排 列 ， 即 自 发 磁 化 至 饱 和 ， 称 为 自 发 磁化 。 铁 磁 体 自 发 磁 化 分 成 若 干 个 小 区 域 （ 这 种 自 发 磁 化 至 饱 和 的 小 区 域 称 为 磁 畴 ） 由 于 各 个 区 域（ 磁 畴 ） 的 磁 化 方 向 各 不 相 同 ， 其 磁 性 彼 此 相 互 抵 消 ， 所 以 大 块 铁 磁 体 对 外 不 显 示 磁 性 。 当温度升高时，原子间距加
4、大，降低了交换作用，同时热运动不断破坏原子磁矩的规则取向，故自发磁化强度下降。单轴各向异性铁磁体具有一个容易磁化的晶轴，原子磁矩的取向只能平行或反平行于这个轴，因此，也只能沿这个轴，或者朝上，或者朝下，由偶然的因素决定。温度升高时，热运动有减弱有序取向的趋势，不过只要温度不太高，仍有为数较多的原子磁矩沿某一取向。直到温度高于居里点，以致完全破坏了原子磁矩的规则取向，自发磁矩就不存在了，材料表现为强顺磁性，其磁化率与温度的关系服从居里外斯定律。为了能更好地体现单轴各向异性铁磁体的顺磁相变过程，我们运用了 Matlab 软件构建了一个简化的伊辛模型，并利用重要性抽样和蒙特卡罗方法模拟了其相变过程。
5、二、建模思想1、伊辛模型（Ising model）Ising 在 1925 年提出了一个描述铁磁体的简单模型，模型虽然十分简单，但是用它讨论铁磁体的相变十分方便。将铁磁体视为 N 个格点组成的 n 维晶格，每个格点均有一个自旋粒子，粒子的自旋计为Si（i=1,2,3,N）它只取+1 或1 两个值，或俗称自旋向上与向下，每个分布描述铁磁系的一个构形。只考虑最近邻自旋相互作用，其作用能的取值原则是：当两个相邻自旋相互平行（沿相同取向）时取-，反平行时，取+；0 对应铁磁性， 0 时，则当所有自旋具有相同取向时系统具有最低的能量，相应于绝对零度下的状态。在足够低的温度下，也会有很多的自旋具有相同的取
6、向，这就是无外场时铁磁具有自发磁化的原因。ii. 物理量的平均值, iillAeiii.物体的磁化率 , 21()MXHkT.物体的磁化强度,ii.物体的热容,221()vcEkT推导过程如下()ssEse热力学统计物理计算机模拟- 4 -2 22()()ss ssEEsEee22()s ss ssEEee,2又:,2kT所以有,22EkT即是,221()vCk四、建立模型下面用 Matlab 语言编写了用 Monte Carlo Method 显示在二维伊辛模型在低温条件下发生相变的简化程序。1、假定有 400 个磁矩，呈 2020 晶格排列，初始状态有 200 个磁矩朝上，200 个磁矩朝
7、下，并在晶格中交替排列，如下图所示（其中 Ising（i,j）=1 表示处于第 i 行第 j 列的晶格点上磁矩朝上并用黑色三角形表示，Ising（i,j）=-1 代表该磁矩朝下，并用红色三角形表示） 。 热力学统计物理计算机模拟- 5 -2、用 rand 语句产生两个随机数，取整并令之等于 i 和 j。改变 Ising（i,j）的符号计算内能 U，并用以下方法判断是否接受：若 ，说明磁 畴 磁 化 方 向 的 改 变 是导致能量体系降低的，在我们的算法里面这样变化我们一定12U接受。若 ，这种情形下，也是接受的。12若 ，说明磁 畴 磁 化 方 向 的 改 变 是 导致体系能量增加的。原则上来
8、说，我们排斥这样的变化，但是系统是存在涨落的，所以尽管能量增加了，这样的改变也会以一定的概率发生。3、随机抽样时要考虑周期性边界条件Ising(1,m)= Ising (21,m)Ising (22,m)= Ising (2,m)Ising (n,1)= Ising (n,21)Ising (n,22)= Ising (n,2)4、本实验中取相对单位,取 K=1，J=1。5、程序流程流程说明：首先建立形如 1 中提到的模型，用 2222 的矩阵表示（实为 2020 的矩阵，四周的元素作为边界条件） ，正负号代表磁化方向，用两个 for-end 语句产生并绘制成图。用 Ising model 的
9、思想建一个二维模型用 monte carlo 方法模拟铁磁-顺磁相变现象先循环 50 万次进行随机翻转再运行 50 万次，统计并计算 Ev,Cv,M,X 的值否则不接受是接受？exp(U1-U2)rand绘制 Ev-T 图,Cv-T 图,M-T 图,X-T 图热力学统计物理计算机模拟- 6 -其次产生两个随机数 sai 和 saj 作为矩阵的行标和列标如 2 所述，将每一组（sai,saj）对应元素翻转（改变正负号） ，再算内能 U 并用 2 中的方法判断是否接受，如此运行 50 万次。接着计算每一次的能量 E及磁化强度 Mo，运行 50 万次，统计得到每个温度对应的 Ev，Cv，M，X，绘制
10、成图。五、结果分析1平均能量与温度的关系图从平均能量与温度的关系图中可以看到随着温度的升高，平均能量 Ev 先逐渐增大最后趋于稳定值，其中在温度（2.0，2.5）间迅速增加，在约为 2.3 处出现拐点，可见相变现象发生在这一温度段中，且转变温度 Tc 约为 2.3。2热容与温度的关系图热力学统计物理计算机模拟- 7 -从热容与温度的关系图中可以看到随着温度的升高，热容 Cv 先增大后减小，在温度约为 2.3 时出现峰值，故转变温度 Tc 约为 2.3。 3磁化强度与温度的关系图从磁化强度与温度的关系图中可以看到随着温度的升高，磁化强度 M 的大小先非常缓慢地减小，大约在温度为 2.35 时迅速
11、减小为零，故转变温度 Tc 约为 2.35。4磁化率与温度的关系图从磁化率与温度的关系图中可以看到随着温度的升高，磁化率 X 先增大后减小，在温度约为 2.4 时出现峰值，故转变温度 Tc 约为 2.4。5. 结果分析综上所述，顺磁物质的平均能量 Ev 随着温度的升高先增大而后趋于稳定值，热容 Cv、磁化率 X 随着温度的升高先增大后减小，磁化强度 M 在转变温度 Tc 处迅速减小为零。这是因为，在温度比较低时，铁磁体内磁畴的有序程度是比较高的，大部分的磁畴都有相同磁化方向，此时系统的内能比较低，但磁化强度较高，随着温度的升高，铁磁体内磁畴的有序程度降低，相应的内能就会增大，而磁化强度则减小，
12、当热力学统计物理计算机模拟- 8 -温度增大到一定值后，变化达到平衡，铁磁体内磁畴的有序程度不再降低，此时内能趋于稳定值，而磁化强度几乎为零。热容是能量对温度的偏导，故在能量图的拐点处热容会出现峰值。综合分析各图，我们得到铁磁-顺磁相变的转变温度约为 2.35。六、讨论总结1. 为使实验模拟结果更有效,更接近实际情况,我们应该注意以下两点：（1）随机取样过程中，选取的样品必须是足够大量的，否则没有统计意义；在本模拟实验中步数取得越大越好。（2）选取的模型必须有代表性，所以随机选取需要有一种非常均匀分布的随机数序列.，所以要求伪随机数发生器要好。2. 做这次课程设计感触还挺大的，理论学起来也许很
13、简单，但是用于实际运用那才是最重要的。这学期的热力学与统计物理和 Matlab 语言的综合运用，让我们感觉到自己所学的知识真的能解决实际问题，当然，这之间经过一个漫长而艰辛的实践历程，多次体会到柳暗花明的感觉。与此同时，我们深深地感到分工合作的重要性，在做课程设计时，我们有的查资料，有的编辑制作，有的编写程序，这么完美的合作才换来一个成功的课程设计。当然，还要感谢老师的细心指导，这让我们少走了不少弯路。七、Matlab 程序clearclc% Ising Model 初始化各磁畴方向Ising=zeros(22,22);for i=1:22for j=1:22if rem(j,2)=0Isin
14、g(i,j)=-1;plot(i,j,rv)hold onelse Ising(i,j)=1;plot(i,j,k)hold onendendendaxis off% Monte Carlo 重要抽样方法使铁磁处于平衡状态Tup=0.1:0.1:5.0;Tdown=5.0:-0.1:0.1;%T 增大for i=1:50E1=0;E2=0;E3=0;M1=0;M2=0;M3=0;热力学统计物理计算机模拟- 9 -for k=1:1000000sai=round(19*rand+2);saj=round(19*rand+2);Ising(1,:)=Ising(21,:);Ising(22,:)=
15、Ising(2,:);Ising(:,1)=Ising(:,21);Ising(:,22)=Ising(:,2);U1=-Ising(sai,saj)*(Ising(sai-1,saj)+Ising(sai+1,saj)+Ising(sai,saj-1)+Ising(sai,saj+1)/Tup(i);U2=-U1;if exp(U1-U2)randIsing(sai,saj)=-Ising(sai,saj);endif k500000 %求出 T 温度下的 Ev,Cv,M,XIsing(1,:)=Ising(21,:);Ising(22,:)=Ising(2,:);Ising(:,1)=Is
16、ing(:,21);Ising(:,22)=Ising(:,2);E=0; Mo=0;for i1=2:21for j1=2:21E=E-Ising(i1,j1)*(Ising(i1-1,j1)+Ising(i1+1,j1)+Ising(i1,j1-1)+Ising(i1,j1+1);Mo=Mo+Ising(i1,j1);endendE=E/2;Mo=Mo/400;E1=E1+E;E2=E2+1; E3=E3+E2;M1=M1+Mo;M2=M2+1;M3=M3+Mo2;endendEv(i)=E1/E2;E2v(i)=E3/E2;Cv(i)=(E2v(i)-Ev(i)2)/(Tup(i)2);
17、M(i)=M1/M2;M2v(i)=M3/M2;X(i)=(M2v(i)-M(i)2)/Tup(i);end% T 减小for i=1:50热力学统计物理计算机模拟- 10 -E11=0;E21=0;E31=0;M11=0;M21=0;M31=0;for k=1:1000000sai=round(19*rand+2);saj=round(19*rand+2);Ising(1,:)=Ising(21,:);Ising(22,:)=Ising(2,:);Ising(:,1)=Ising(:,21);Ising(:,22)=Ising(:,2);U11=-Ising(sai,saj)*(Ising(
18、sai-1,saj)+Ising(sai+1,saj)+Ising(sai,saj-1)+Ising(sai,saj+1)/Tdown(i);U21=-U11;if exp(U11-U21)randIsing(sai,saj)=-Ising(sai,saj);endif k500000 %求出 T 温度下的 Ev,Cv,M,XIsing(1,:)=Ising(21,:);Ising(22,:)=Ising(2,:);Ising(:,1)=Ising(:,21);Ising(:,22)=Ising(:,2);E0=0;Moo=0;for i1=2:21for j1=2:21E0=E0-Ising(i1,j1)*(Ising(i1-1,j1)+Ising(i1+1,j1)+Ising(i1,j1-1)+Ising(i1,j1+1);Moo=Moo+Ising(i1,j1);endendE0=E0/2;Moo=Moo/400;E11=E11+E0;E21=E21+1; E31=E31+E02;M11=M11+Moo;M21=M21+1; M31=M31+Moo2;endendEv0(i)=E11/E21;E2v0(i)=E31/E21;Cv0(i)=(E2v0(i)-Ev0(i)2)/(Tdown(i)2);
