铁电体相变三维Ising模型的Monte Carlo 模拟

一格点 Si 只考虑最近邻的六个格点上的自旋 S1_ S6 ,且仅考虑单轴各向异性情
形，即所有的 Si 只能取±1 。则在无外场情况下，上述伊辛模型可简化为：

H
=
−
1 2
∑
i, j
J ij
Si
S
j
（）
exp(-6J/kT)
4J
exp(-4J/kT)
2J
exp(-2J/kT)
0
1
2J
1
4J
1
6J
1
C.计算机模拟过程
（1） 选定点阵规格 /×/×/，由于格点数越多，所需时间越长，仅对 /，，
三种情况进行讨论；取 + L5
的起始值。 （2） 任取一个自旋点阵排列作为起始状态，由每一个格点 SJ 及其近邻自旋情
关键词：Ising 模型 Monte Carlo 模拟 居里温度 极化强度 自旋组态
Monte－Carlo Study of Three Dimensional Ising Model
in Ferroelectric Phase Transition
Abstract: The spontaneous polarization and phase transition in order-disorder ferroelectric crystal was studied using Monte Carlo computer simulation. Our study has been developed in the frame work of the three dimensional Ising model of From our results, it is concluded that, in agreement with the analysis of statistical mechanics, the polarization changed significantly near the ferroelectric phase transition Curie temperature and the Curie temperature was different because of the different scale. The spin configuration of this model at various temperatures was also obtained. Key words: Ising model; Monte Carlo simulation; Curie temperature; polarization; spin configuration
H
=
∑
i
H (i) 1
+
1 ∑
2 i, j
H2 (i,
j)
(1)
又因为系统的运动可以用粒子的泡利自旋算符 S x, S y , S z 来表示。所以哈密顿
量亦可以写成【1】：
∑ ∑ H = −Ω i
Si x
−
1 2
i, j
J ij Si z S j z

(2)
图 1 双势阱模型
图 2 三维自旋格点模型
P
=
⎪⎨⎧exp(−
∆H kT
)
∆H ≥ 0
⎪⎩ 1
∆H < 0
（5）
对所有格点均进行这样的处理后，系统的状态即朝着某温度 T 对应的平衡态前进
了一步，称之为一个蒙特卡罗步（MCS）.经过足够多的蒙特卡罗步，系统就进入
平衡态。上述过程实际上就是构造系统状态演化的马尔可夫（Markov）过程，因
此选取的 SJ 反转概率必须满足细致平衡条件的要求。表 1 给出了三维情况时 SJ

由表中数据不难看出，随着温度的升高，极化强度亦呈上升的趋势<>，而且
在某一温度区间（即 + L5
_ ）出现了大幅增加的情况<>，由铁电
相变理论知，居里温度就在该区间内<>。如果希望更为准确地确定居里点的位置，
则可通过缩小 + L5
的取值范围和增加其取值密度（减少其增量）再进行模拟计
况，按表 1 计算出 SJ 反转的概率 P。
（3） 若 Q 则接收 SJ 的反转；若 Q ≤ 则由<，> 产生一个均匀分布随机数
ξ，与得到的 P 值比较：若 P>ξ则接受 SJ 的反转；若 P≤ξ则认为 SJ 保 持不变。 （4） 按（2）和（3）确定点阵上所有格点的自旋指向，即完成了一个蒙特卡罗 步（.$4）。 演化足够多的 MCS 后，系统达到平衡态。选取点阵格点数越多，需要的 .$4 也越多。每隔若干 MCS 抽取一个状态作为样本，计算样本的每格点极化强度 Ps = ∑ Si N 3 ,求和对所有格点进行[5]。收集足够多的样本，求平均每格点极化强
算。用这种方法计算居里温度过程简单，效率高，可以获得较为精确的数据。
同时，由表中的数据还可以看出，对三种不同规格的体系来说，极化强度发生突 变的温度范围不一样。理论上，随着体系尺寸的增加，其居里温度 5D 呈单调上 升趋势。表 中的数据与之基本相符，但并不能明显看出规格 和规格 之间 5D 的差异，如上述，若进一步缩小 + L5
自旋的四体相互作用，从而统一说明了二级相变和一级相变，并使自发极化强度 与温度关系的计算结果与实验结果符合得很好；曲保东等采用平均场近似研究了 铁电超晶格、铁电体薄膜的相变问题。伊辛模型二维情况的严格解在 1944 年由 昂萨格求出，但三维情况的严格求解仍是当前富有挑战性的研究课题。而 .POUF $BSMP 方法提供了一种非常有效的途径。因此，采用这种方法运用计算机模拟， 通过分析体系自发极化强度随温度的涨落行为，能较准确地找到铁电体系居里点 的位置，并且可以直观地显示 *TJOH 自旋点阵指向分布随温度的下降而由无规向 有序的演化过程。
1 引言
铁电体铁电相与顺电相之间的转变通常简称为铁电相变，转变的临界点温度 称为居里温度或居里点 Tc。铁电相变的实质是出现自发极化，且自发极化可以 在电场作用下发生转向。根据晶体结构测定和理论分析，可将铁电相变分为位移 型和有序无序型。其中，对有序无序型铁电相变，其微观理论主要是赝自旋理论， 而横场 *TJOH 模型即该理论中的一个重要模型。很多方法在近似处理该模型体系 是相当成功的：#MJOD 首次采用平均场近似计算了铁电体的赝自旋平均值，成功 地解释了顺电铁电相变及其临界温度的同位素效应；此后 (FOOFT 等提出了描述 有序无序型铁电相变的横场伊辛模型；王春雷等用格林函数理论成功的求解了赝
求和只对每个格点的最近邻格点进行。对于点阵边缘的格点，采用周期性的边界
条件确定 SJ 的近邻。即：如果 SJ 在立方体上底面，则取立方体的下底面为其上
各个格点的相应上近邻，最左（后）面格点则取最右（前）面相应格点为其左（后）
近邻，如此等等。
在任一温度下，不存在外电场时，SJ 的指向受两种因素影响。一方面，自旋

规格 2
规格 3
间的相互作用使近邻自旋指向趋于规则排列，以保持能量最低。规则排列的自旋 使极化强度不为零，这是自发极化的根本原因。由 SJ 及其近邻格点上的自旋指 向，计算 SJ 反转前后的哈密顿量。反转前，H1= -JÔSJSK -JSJ S S … S
反转后，H2= -JÔSJSK -JSJ S S … S
分别是第
J
个质子自旋算符的
Y
分量，[
分量和第
K
个质子自
旋算符的 [ 分量，+JK 是相互作用系数，μ 为有效偶极矩。若将 Ω 视为横向场，
J 视为交换积分，则（3）式即横场伊辛模型。对该模型进行适当简化之后，进
行 .POUF$BSMP 计算机模拟。
3 .POUF$BSMP 模拟
具体做法是： " 模型的简化 假设三维自旋网格是一个 N×N×N 的立方体点阵，对网格中的任
所有 7 种近邻情况及其相应的反转概率，其中，哈密顿量由（4）式算出，而概
率 P 由（5）式确定。x 表示自旋指向与 SJ 相同的最近邻格点数。
表 1 三维格点自旋的近邻情形及相应反转概率（取演化前 SJ）
X
H1
H2
6
-3J
3J
5-2J2J4-JJ
3
0
0
2
J
-J
1
2J
-2J
0
3J
-3J
ΔH
P
6J
其中 S z 量度了左右平衡位置上粒子占据数之差亦即量度了有序化的程度。Ω 为
隧穿频率，J 是耦合常数。当体系加一外场 E 时，
H
= −Ω∑ Si x i
−
1 2
∑
i, j
J
ij
Si
z
S
j
z
− 2µ ∑ ESiz i
（）
式中，
Six
,
S
z i
,
S
z j
。其中 SJ 、SJ'分别代表反转 前后格点 SJ 的自旋 设 SJ =1,则 SJ' = -1。当 H2 > H1 时，格点从自旋 SJ 转向 到 SJ'比较困难，当 H2 < H1 时，反转容易发生。另一方面，热运动效应则使自 旋指向趋于无规，当温度较高时，即使 H2 > H1 反转也有一定的概率发生。两 种因素的竞争导致临界点的存在：当温度小于某个值 5D 时，自旋点阵中各格点 的自旋指向呈现有序；一旦超过 5D 则处于无规状态，极化强度为零。 # 蒙特卡罗算法 从点阵系统中随机的产生一初始状态，并假设特定温度 T。 从该初态开始，按一定抽样方法决定每个格点的自旋是否反转。设Δ)H2 - H1 表示 SJ 反转前后哈密顿量的改变量，SJ 发生变化的概率 P 满足[2]：
i
度〈Ps〉。 （5） 将 + L5
取一个新的值，重复步骤（）_（），即可得到〈Ps〉随 + L5
值的变化关系。
结果与分析
（）下面的表 给出了 ××，××，×× 点阵 用规格 _ 表
示
的计算结果。
表 三种规格的三维点阵的<Ps>随 + L5
2 *TJOH（伊辛）模型
伊辛模型是相变统计理论的一个核心模型，它是建立在双势阱模型（如图 1） 基础之上的。在有序无序型铁电晶体中，导致铁电相变的离子运动是双势阱间的 运动。这类晶体中含有氢键，顺电相时氢（质子）在氢键中两个可能位置上等概 率分布，呈无序状态，而在铁电相时氢择优地占据这两个可能位置之一，呈有序 状态。因此，氢的有序化是有序无序型铁电相变的触发机制。如图示，质子处在 两个势阱之中，在一定条件下可以贯穿势垒，由一个阱进入到另一个阱。借用铁 磁理论中的自旋波理论，设想每一个这样的单元用一个赝自旋代表。质子位于左 右两个势阱相应于赝自旋的上下两种取向，整个晶体中质子的分布和运动则用系 统的赝自旋波来描述。在计入贯穿势垒的隧道效应、忽略高能级的状态及粒子在 阱内的运动的情况下，赝自旋系统的哈密顿量显然包括单粒子部分 H1(i)和相互 作用部分 H2(i, j),即
《计算材料学》课程设计
指导老师：江建军 教授 电子科学与技术系 2004 年 6 月
铁电体相变三维 Ising 模型的 Monte Carlo 模拟
蔡梦 曹亨 何祖辉 蓝天 李家岗 林键 孙文建 余会星 余志龙 张鲁 张君
摘 要:本文采用横场 Ising 模型和 Monte Carlo 方法对有序无序型铁电体的自发 极化和其相变特征进行了计算机模拟。研究结果表明，在铁电体发生相变的临界 点温度即居里温度附近，极化强度发生强烈的改变，且居里温度随模型尺寸的大 小不同而异。还得到了该模型不同温度下自旋组态图样，完全符合统计力学的分 析。
的变化
+,5 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.40 0.41 0.42 0.43
规格 1