（凝聚态物理专业论文）ising系统零温动力学的monte+carlo研究.pdf

摘要 本文采用m o n t ec a r l o 方法研究了二维i s i n g 系统从高温淬火至零温的g l a u b e r 动力 学演化过程，主要内容包括以下几个方面： 1 研究了系统初始位形的状态对系统演化的影响。对于初始磁化强度为零的位 形( 即方向向上和向下的自旋个数相等) ，发现演化终态( 基态或条形态) 具有某种对称 性，如终态的磁化强度关于零对称；自旋方向向上的基态、方向向下的基态和条形态所 占比例相同。对初始位形中磁化强度不等于零( 即不同方向的自旋个数不相等) 的情况， 这种对称性会消失。 2 研究了系统长时演化出现的反常行为。在演化到一定阶段后，系统会出现一个极 长的演化过程，导致一些物理量出现反常行为。如能量演化和”持续性问题所遵守的幂 律衰减会被打破，而出现一个演化过程中的平缓期。通过研究反常行为中磁化强度和能 量的变化，发现磁化强度在整个演化过程中一直在波动，不能作为表征该反常行为的物 理量，而系统能量在反常过程中保持不变，可以用此稳定能量值来表征系统演化中的反 常行为。通过对比不同大小的晶格在反常行为发生时所对应的系统能量，发现对于基态 实现( r e a l i z a t i o n ) ，该能量值为4 l ，而对于条形态实现，该能量值接近2 l ，这里l 是 晶格大小。 3 从基态和条形态位形入手，分析了动力学反常过程的原因。在基态实现中，对角 线型位形结构的特殊性是造成基态演化出现反常行为的原因，而在条形态实现中，其极 长的演化阶段发生在演化末期，所出现的锯齿状位形是造成条形态实现反常的原因。计 算了这两种特殊位形的能量，所得结果与模拟结果相符。 4 计算了系统终态的能量。二维i s i n g 系统终态只有两种一基态和条形态，基态能 量为o ，条形态中磁畴分界线条数绝大多数是2 ，其能量约等于2 l 。根据基态和条形态 的能量以及两者在模拟中出现的比例，系统终态的能量近似等于2 3 l 。 5 研究了m e t r o p o l i s 算法和c o n s t r a i n t 算法对系统演化的影响。和g l a u b e r 动力学 比较，发现m e t r o p o l i s 算法对动力学过程影响不大，而c o n s t r a i n t 算法会使动力学结果 产生明显变化。 关键词：i s i n g 系统；m o n t ec a r l o 方法；g l a u b e r 动力学；零温动力学 a b s t r a c t u s i n gm o n t ec a r l om e t h o d ，w ei n v e s t i g a t et h ed y n a m i c a le v o l u t i o no ft w o d i m e n s i o n a l i s i n gs y s t e mq u e n c h e df r o mh i g ht e m p e r a t u r et oz e r ot e m p e r a t u r e t h em a i nr e s u l t so ft h i s t h e s i sa r ea sf o l l o w s ： 1 t h ee f f e c to fi n i t i a lc o n f i g u r a t i o no ft h ei s i n gs y s t e mo ni t sd y n a m i c a le v o l u t i o ni s i n v e s t i g a t e d i ft h em a g n e t i z a t i o ni nt h ei n i t i a lc o n f i g u r a t i o ne q u a l st oz e r o ( i nt h i sc a s et h e n u m b e ro fu ps p i n si se q u a lt ot h a to fd o w ns p i n s ) ，t h ef i n a ls t a t ei ss y m m e t r i c f o re x a m p l e ， t h em a g n e t i z a t i o ni ss y m m e t r i ca b o u tz e r oa n dt h er a t i oo f u p s p i ng r o u n ds t a t ei st h es a m ew i t h d o w n - s p i ng r o u n ds t a t eo rs t r i p es t a t e i ft h em a g n e t i z a t i o ni nt h ei n i t i a lc o n f i g u r a t i o nn o te q u a l t oz e r o ，t h es y m m e t r yw i l ld i s a p p e a r 2 t h ea n o m a l o u sb e h a v i o ri nl o n g t i m ed y n a m i c a le v o l u t i o ni ss t u d i e d a tt h el o n g t i m e s t a g eo ft h ed y n a m i c a le v o l u t i o n ，t h e r ei sav e r yl o n gr e l a x a t i o nw h i c hc a u s i n gs o m ea n o m a l o u s b e h a v i o r s f o re x a m p l e ，t h ep o w e r - l a wd e c a yi ne n e r g ya n dp e r s i s t e n c ew i l lb eb r o k e n b y s t u d y i n g t h ef l u c t u a t i o n so fm a g n e t i z a t i o na n de n e r g yi nt h ea n o m a l o u sd y n a m i c s ，w ef i n dt h a t t h ee n e r g yi ss t e a d ya n dc a l lb eu s e dt oc h a r a c t e rt h ea n o m a l o u sb e h a v i o r , b u tt h em a g n e t i z a t i o n a l w a y sf l u c t u a t e s ，t h u si tc a nn o tb eu s e dt oc h a r a c t e rt h ea n o m a l o u sb e h a v i o r i na d d i t i o n ，w e f i n dt h a tw h e na n o m a l o u sb e h a v i o ra p p e a r st h ee n e r g yi s4 la n ds l i g h t l yl a r g e rt h a n2 li n + g r o u n ds t a t er e a l i z a t i o na n dt h es t r i p es t a t er e a l i z a t i o n ，r e s p e c t i v e l y , w h e r el i st h el a t t i c es i z e 3 b ys t u d y i n gt h ec o n s t r u c t i o n so ft h eg r o u n ds t a t ea n ds t r i p es t a t e ，w ef i n dt h es o u r c e s o fa n o m a l o u sd y n a m i c a le v o l u t i o n i ng r o u n ds t a t er e a l i z a t i o n ，t h ea n o m a l yo r i g i n a t e df r o m t h ed i a g o n a ls t r i p ec o n f i g u r a t i o n ，a n di ns t r i p es t a t er e a l i z a t i o n ，t h ea n o m a l yc o m e sf r o mt h e s a w t o o t h l i k ec o n f i g u r a t i o nw h i c ha p p e a r si nt h el a s ts t a g eo fe v o l u t i o n t h ee n e r g yv a l u e so f t h es p e c i a lc o n f i g u r a t i o n sa r ei ng o o da c c o r d a n c ew i t ht h o s ei ns i m u l a t i o n 4 t h ee n e r g yo fs y s t e mi nf i n a ls t a t ei sc a l c u l a t e d i ti sf o u n dt h a tf o r t h et w o d i m e n s i o n a l i s i n gs y s t e m ，t h e r ea r et w of i n a ls t a t e s - - g r o u n ds t a t ea n ds t r i p es t a t e t h ee n e r g yo fg r o u n d s t a t e i s0 ，a n dt h a to fs t r i p es t a t ei sa p p r o x i m a t e l y2 lb e c a u s et h en u m b e ro fb o u n d a r i e sb e t w e e n d o m a i n sm a i n l yi s2 i ti sa l s of o u n dt h a tt h ee n e r g yo ft h es y s t e mi nf i n a ls t a t ec a nb ew r i t t e n a s2 3 l 5 t h ee f f e c to ft h em e t r o p o l i sa l g o r i t h ma n dc o n s t r a i n to n eo nd y n a m i c a le v o l u t i o ni s a l s od i s c u s s e d c o m p a r i n gw i t hg l a u b e ra l g o r i t h m ，w ef i n dt h a tt h ee f f e c to ft h ec o n s t r a i n t a l g o r i t h mi ss t r o n g e rt h a nt h a to ft h em e t r o p o l i sa l g o r i t h m k e yw o r d s ：i s i n gs y s t e m ；m o n t ec a r l om e t h o d ；g l a u b e rd y n a m i c s ；z e r o t e m p e r a t u r e d y n a m i c s 曲阜师范大学硕士学位论文原创性说明 ( 在口划“) 本人郑重声明：此处所提交的博士口硕士面论文i s i n g 系统零温动力学的m o n t e c a r l o 研究，是本人在导师指导下，在曲阜师范大学攻读博士口硕士留学位期间独立 进行研究工作所取得的成果。论文中除注明部分外不包含他人已经发表或撰写的研究成 果。对本文的研究工作做出重要贡献的个人和集体，均已在文中已明确的方式注明。本 声明的法律结果将完全由本人承担。 作者躲戎描 日期：彩，石f o 曲阜师范大学硕士学位论文使用授权书 ( 在口划“，) i s i n g 系统零温动力学的m o n t ec a r l o 研究系本人在曲阜师范大学攻读博士口 硕士囱学位期间，在导师指导下完成的博士口硕士囱学位论文。本论文的研究成果归 曲阜师范大学所有，本论文的研究内容不得以其他单位的名义发表。本人完全了解曲阜 师范大学关于保存、使用学位论文的规定，同意学校保留并向有关部门送交论文的复印 件和电子版本，允许论文被查阅和借阅。本人授权曲阜师范大学，可以采用影印或其他 复制手段保存论文，可以公开发表论文的全部或部分内容。 作者签名：戒枣狰 新撕2 撕友d 代 日期：如矽f f o 黾强：加g 、乇、k 第一章绪论 人们早在三千多年前就发现一些物体具有磁性，但磁学理论的研究在十九世纪初才 得到迅速发展。o e r s t e d 在1 8 2 0 年的发现，使人们认识到电子的运动一电流可以产生磁 场。1 9 2 5 年，g e u h l e n b e c k 和s g o u d s m i t 假设电子具有自旋【1 】，并且其行为类似于 小磁棒，在外磁场中电子磁场的方向与外磁场方向相同或相反。同一年，l e n z 建议他的 学生e i s i n g 考虑自旋间的相互作用，则晶格上平行的自旋相互吸引，反平行的自旋相 互排斥，这就是i s i n g 模型的由来。自从e i s i n g 于1 9 2 5 年提出i s i n g 模型以来 2 】，i s i n g 模型已成为统计物理学中重要的理论模型之一。m o n t ec a r l o 方法作为常用的模拟方法 之一，在经典粒子系统、自旋系统、分形系统以及逾渗、聚合物、晶体生长等物理问题中 得到了广泛应用 3 ，4 】。本文采用m o n t ec a r l o 方法对i s i n g 自旋系统的零温动力学进行 研究。 1 1 自旋系统与i s i n g 模型 研究自旋系统首先要确定系统位于什么样的晶格结构上，以此来决定自旋之间相互 作用的形式和大小 5 】。晶格结构一般由线和顶点( v e r t e x ) 构成。晶格结构中的顶点，称 为座点( s i t e ) 。座点上可以放置自旋，也可以空置。座点之间的连线称为晶格键( 1 a t t i c e b o n d ) 。 图1 1 方晶格结构 本文所研究的晶格是方晶格，是最常见的晶格结构之一，如图1 1 。在大部分研究工 作中，涉及到的晶格结构都是指所有相邻座点之间的晶格键都存在，并且所有的座点都 被自旋占据，这样的晶格称为纯晶格。还有些研究工作中的晶格结构其相邻座点间只存 在部分晶格键，这种情况称为键稀释晶格( b o n d - d i l u t e dl a t t i c e ) ；或只有部分座点被自旋 占据的情况，称为座点稀释晶格( s i t e d i l u t e dl a t t i c e ) 。 第一章绪论 i s i n g 系统作为最简单的自旋系统之一，其哈密顿量可以表示为 h = 一j t j a i a j ， ( 1 1 ) ( 巧) 其中c r i 代表格点i 处的自旋，自旋取值+ 1 或一1 ，五f 是格点i 和j 上的自旋相互作用 耦合系数，求和遍历最近邻自旋对。对铁磁模型，五， 0 ，反铁磁模型也， 0 。每个自 旋上的磁化强度定义为 m = 寺f 吼， ( 1 2 ) m 2 丙乙吼， 【l z ) 。 t 其中是自旋总数，单个自旋磁化强度的值在+ 1 和一1 之间。 e i s i n g 只计算了简单一维链上的i s i n g 模型，发现只有温度t = 0 时才会发生相 变【2 】。l a r so n s a g e r 研究了方晶格上的二维i s i n g 模型 6 ，得到了严格解并发现存在相 交，该结果被视为临界现象发展史上的一个里程碑。由于运算的复杂性，三维i s i n g 模型 的严格解一直没有得到。现在相变理论的研究，多采用各种近似方法，如l a n d a u 平均场 理论，重整化群理论，m o n t ec a r l o 方法等。 1 2 自旋系统的动力学 根据所研究自旋系统所处的温度不同，动力学常分为临界动力学和零温动力学。临 界动力学主要研究系统在临界点附近发生的非平衡现象，即研究临界点附近系统由非 平衡态向平衡态的过渡过程。在演化过程中，系统会出现一些复杂而有趣的现象，如描 写动力学性质的物理量在临界点附近会发散，临界动力学就是从理论上来分析这些现 象【7 】。零温动力学主要研究系统从高温淬火至零温时系统在演化中所表现出来的规律。 m o n t ec a r l o 方法是研究动力学理论的一种重要方法，所采用的算法主要有单自旋 翻转、集团翻转等。在单自旋算法中，常见的有两种：单个自旋反转和最近邻自旋交 换【8 】。单个自旋反转算法没有守恒律，会改变系统的总磁化强度，实现的是一个巨正则 系综，又称为g l a u b e r 动力学 9 】；最近邻自旋交换算法中总磁化强度守恒，实现的是一 个正则系综，又称为k a w a s a k i 动力学 1 0 。 1 3 课题研究历史及现状 1 9 9 4 年，d e r r i d a 等人研究了i s i n g 模型和p o t t s 模型在零温动力学中的持续( p e r - s i s t e n c e ) 现象，该现象是指从开始模拟直到时间t 一直没有翻转的自旋占总自旋的比 例【11 。他们发现持续现象遵循幂律衰减 p ( t ) 一t ， ( 1 3 ) 其中p ( t ) 代表t 时刻未翻转自旋的比例。对二维i s i n g 模型来说，持续指数0 0 3 7 ，并 认为该指数具有非平凡性( n o n t r i v i a l ) 。随后，s t a u f f e r 在更大系统上计算了一维至五维 2 第一章绪论 i s i n g 模型的持续指数【1 2 】，验证了d e r r i d a 等人的结果。1 9 9 5 年，d e r r i d a 将一维p o t t s 模 型类比反应一扩散模型( r e a c t i o n d i f f u s i o nm o d e l ) ，得到持续指数的解析解 1 3 1 。随后，持 续问题得到了广泛研究 1 4 - 1 8 】。 在零温动力学中，另一个感兴趣的问题就是i s i n g 铁磁系统演化的最终状态。对于 无限大系统，随着系统的演化，磁畴以亡1 2 的速度增长【1 9 ，自旋可以进行无限次翻 转【2 0 。而对于m o n t ec a r l o 模拟而言，采用的系统都是有限大小，有限尺寸效应不容忽 视。2 0 0 1 年，s p i r i n 等人的模拟结果发现系统最终会冻结在某种状态：基态( g r o u n ds t a t e ) 或亚稳态( m e t a s t a b l es t a t e ) ( 见图1 2 ) 。基态是指所有的自旋方向都一致，都向上或者都 向下。处于亚稳态的二维i s i n g 模型也不再进一步演化，其位形由宽度大于2 的交替水 平( 竖直) 自旋条构成，也称为条形态。终态中基态和条形态所占的比例是2 ：1 ，并且与 模拟采用的晶格大小无关【2 1 ，2 2 1 。同一年，s v e n s o n 研究了i s i n g 模型在随机图上的冻 结现象【2 3 1 。2 0 0 5 年，s u n d a r a m u r t h y 和s t e i n 进一步研究了初始位形对终态的影响 2 4 。 设初始位形中自旋值取+ 1 的概率是p ，模拟结果显示当p = 1 2 时，到达条形态的比例 达到极大值；当p 1 2 时，绝大部分初始位形演化到基态，而只有一少部分演化到条形 态，p 偏离1 2 越大越明显。对于p 1 2 情况，结果还显示对无限大系统，系统的演化 也是有限的，即冻结到某个终态而不再演化。 + + 图1 2 基态和条形态示意图：( a ) 基态；( b ) 条形态。+ 代表白旋方向向上，一代表自旋方向 向下 在零温动力学研究中，某些位形具有反常的演化过程。1 9 9 9 年，l i p o w s k i 发现线性 长度为三的磁畴的平均寿命与l 2 成正比，但有些初始位形的平均寿命与l 3 成正比，他 认为这种长寿命位形的存在是造成相有序( p h a s eo r d e r i n g ) 动力学反常的原因【2 5 1 。2 0 0 1 年，s p i r i n 等人认为零温动力学驰豫过程由两个截然不同的时间尺度构成。对于方晶格 上的二维i s i n g 模型，其长时间的驰豫过程来源于长寿命的对角线型条形( d i a g o n a ls t r i p e ) 缺陷 2 1 1 。2 0 0 5 年，d eo l i v e i r a 等人根据终态的不同，分别统计了基态实现( r e a l i z a t i o n ) 和条形态实现，进一步证实这两种实现由两种时间尺度主导 1 5 。 3 + + + + + + + + + + + + + 十 + + 第一章绪论 1 4 本论文研究的问题及意义 前面已经提到，零温动力学中二维i s i n g 模型的终态有两种：基态和条形态，并且这 两种状态实现过程中具有不同的驰豫时间尺度。在本文中，我们从基本的位形演化入手， 采用m o n t ec a r l o 方法进一步研究i s i n g 模型动力学演化中出现的这两种实现过程，主要 内容如下： 1 研究了造成极长驰豫时间的原因。该过程的了解，有利于在模拟中提高效率，改 进算法，并为研究反常动力学演化提供依据。 2 特殊位形的能量表达式。我们通过观察系统在演化中的能量变化，可以判断系 统演化所处的阶段，从而对整个演化过程得以更清楚的了解。 3 终态位形的能量。终态能量在磁畴的计算中占据重要的地位，因此该能量表达 式对于研究磁畴相关问题具有重要意义。 4 第二章m o n t ec a r l o 方法 m o n t ec a r l o 方法和分子动力学方法( m o l e c u l a rd y n a m i c sm e t h o d ) 是数值模拟中两种 主要的计算方法【2 6 】。分子动力学方法则是根据经典力学从原子、分子水平上求解多体 问题，但运算量大，程序复杂。m o n t ec a r l o 方法是一种随机统计方法，对于与外界存在 复杂相互作用的系统或系统本身就很复杂的问题，此时求解微观运动方程比较困难，而 m o n t ec a r l o 方法可以给出近似解，因此得到广泛应用。 2 1随机数 在m o n t ec a r l o 方法中，随机数问题占据着重要的地位。真正随机的数字序列不容 易得到，而且不容易重复。因此，现在的随机数一般采用代数运算的方法获得，得到的 随机数具有确定性，称为伪随机数( p s e u d o - r a n d o mn u m b e r ) ，下面提到的随机数均指伪随 机数。随机数质量的好坏直接影响到模拟结果的精度，随机序列之间的关联会直接导致 模拟出现错误的结果。比如，1 9 9 2 年f e r r e n b e r g 等人发现：一些被认为是高品质的随机 数在w o l f f 集团翻转算法中得出了错误的结果【2 7 。 最早期的随机数是由j o h nv o nn e u m a n n 提出的“中间_ 平方”法得到的，该方法曾用 于原子弹的研究之中。中间_ 平方法的运算公式为k + 1 = m i d d l e _ d i g i t s ( 墨k ) ，x 代 表随机数。假设我们取六位的随机数，前一个随机数为6 5 4 3 2 1 ，将它平方后得到 4 2 8 1 3 5 9 7 1 0 4 1 ，取中间六位，则新的随机数是1 3 5 9 7 1 。这种方法产生的数字序列，虽然 看起来是随机的，其实并不随机。现在常用的随机数发生器有线性同余法、移位寄存器 序列以及不同方法的组合等【2 8 3 0 。 ( 1 ) 线性同余法 线性同余法是由l e h m e r1 9 5 1 年提出的，由于计算简单，所以得到了广泛应用。其 递推公式为一 k + 1 = ( a x + c ) m o dm ， ( 2 1 ) 其中a 、c 、m 分别称为乘子、增量和模，m o d 代表程序中的模运算。这些参数的值并不 能任意选择，否则会产生不随机的序列。比如，我们取m = 1 0 ，) c o = a = c = 7 ，则所产 生的序列为 7 ，6 ，9 ，0 ，7 ，6 ，9 ，0 ， 是一个循环序列。因此参数的选择应该遵循一定的规则，以保证随机序列的质量。该方 法中的参数如果取 + 1 = 1 6 0 8 7 m o d ( 2 3 1 1 ) ， ( 2 2 ) 5 第二章m o n t ec a r l o 方法 就是著名的1 6 8 0 7 随机数发生器，周期可以达到2 3 1 2 。另外一个常用的选择是 + 1 = 6 9 0 6 9 m o d2 3 2 ，( 2 3 ) 周期可达到2 3 0 。线性同余发生器的缺点是周期不够长，并且具有明显的高维相关性。 ( 2 ) 移位寄存器序列 移位寄存器序列发生器又称t a u s w o r t h e 序列发生器，是t a u s w o r t h e 于1 9 6 5 年提出 来的。它基于多项式与模2 运算，又称异或运算，符号记作o 。例如多项式为x p + x q + 1 ， 则有对应的移位寄存器序列发生器为 = ( k 呻+ k 一( p 一口) ) m o d2 ，竹= p ，p + 1 ， ( 2 4 ) 或者 = k _ p0 七一口) ，n = p ，p + 1 ， ( 2 5 ) 著名的随机数发生器“r 2 5 0 ”采用p = 2 5 0 ，q = 1 0 3 ，共需要2 5 0 个初始的随机数。 前面提到，线性同余法具有高维相关性，因此有人提出了非线性同余法。非线性同 余法在程序实现中使用的i f 判断语句，不能在超级计算机上进行向量化并且不适合高速 m o n t ec a r l o 模拟，因此移位寄存器法在提出时得到大量应用。但不少人认为这是个坏的 随机数发生器。前面提到f e r r e n b e r g 等人在对i s i n g 模型研究采用集团翻转算法时发现 了严重的错误，其采用的随机数发生器就是r 2 5 0 。另外，g r a s s b e r g e r 于1 9 9 3 年在三维 自回避行走中也发现t a u s w o r t h e 序列发生器是一种不好的随机数发生器【3 1 。 ( 3 ) k i s s 算法 k i s s ( k e e p i ts i m p l e ，s t u p i d ) 随机数发生器是m a r s a g l i a 提出来的一种组合发生器。 它组合了三个简单发生器：同余( c o n g m e n t i a l ) 发生器、三移位( 3 - s h i f t ) 发生器和进位乘 ( m u l t i p l y w i t h c a r r y ) 发生器。其c 语言源代码如下： u n s i g n e dl o n gk i s s ( v o i d ) s t a t i cu n s i g n e dl o n gx = 1 2 3 4 5 6 7 8 9 ，y = 3 6 2 4 3 6 ， z = 5 2 1 2 8 8 6 2 9 ，c = 7 6 5 4 3 2 1 ； u n s i g n e dl o n g l o n gt ，a = 6 9 8 7 6 9 0 6 9 l l ； x = 6 9 0 6 9 * x + 1 2 3 4 5 ； y “= ( y 1 7 ) ；y “= ( y 3 2 ) ； r e t u r nx + y + ( z = t ) ； h t t p ：o l d m i l l u c h i c a g o e d u w i l d e r c o d e r a n d o m p a p e r s m a r s a g l i a _ 2 0 0 3 h t m l 6 第二章m o n t ec a r l o 方法 由于k i s s 发生器简单，速度快，通过了所有随机数测试，且周期达到2 1 2 4 ，是被广 泛应用的随机数发生器之一。 ( 4 ) 马其赛特旋转( m e r s e n n et w i s t e r ) 马其赛特旋转是m a t s u m o t o 和n i s h i m u r a 于19 9 8 年提出的【3 2 ，也是我们模拟中所 采用的随机数发生器。该方法依据如下线性递归 x k + n = x k + mo ( j 罐i ) 罐+ 1 ) a ，七= 0 ，1 ，( 2 6 ) 其中参数的含义请参考文献【3 2 。 该发生器具有2 1 9 9 3 7 1 的超长周期，并且没有高维相关性；在模拟过程中运算速度 快，所需时间短，并通过了所有的随机数测试，被认为是最好的随机数发生器之一。 另外还有一些其他的随机数发生器，如经典和延迟f i b o n a c c i 序列等。虽然不同的随 机数发生器有着这样或那样的缺点，以至于有人认为“寻找一个真正随机的伪随机数发 生器是一个不可能完成的任务”【3 3 ，但对本文的模拟而言马其赛特旋转随机数发生器 足以满足我们的要求。 2 2 m o n t ec a r l o 方法简述 2 2 1引言 我们先举例来说明m o n t ec a r l o 方法的中心思想。比如计算常数7 r 的值，比较著名 的有b u f f o n 投针试验，这里我们采用另一种方法：首先计算1 4 圆的面积s ( 见图2 1 ) ， 再根据公式4 s n 2 求丌，r 是圆的半径。为方便，我们取r = 1 。首先我们连续生成范 围在 0 ，1 】之间的随机数对，对应于z 和可轴。每一对随机数对应一个坐标点( z ，可) ，计 图2 1 根据圆面积求7 r 值 7 第二章m o n t ec a r l o 方法 算? - 2 = z 2 + y 2 ，如果7 t 2 冬r 2 ，将计数器m ( 用来统计落入圆内点的个数) 加1 。当产生 个坐标点后，落入1 4 圆内部的点个数与总个数之比m n 近似等于圆的面积s ，则 7 r = 4 m n r 2 ，( 2 7 ) 从上面可以看出，m o n t ec a r l o 方法的中心思想就是首先根据所求问题构造合适的 模型，让待求物理量成为该模型的一个参数，然后利用随机变量的数值序列求得该问题 的近似解。该方法最早用于具有随机性的中子散射和吸收过程的研究，目前已广泛用于 数学积分，统计物理和量子力学等学科的计算之中。 s = 击、r 2 一z 扎 ( 2 9 ) ，= ( 6 一。) 击m 1 ) - ( b - a ) ( 玑 ( 2 1 1 ) ，兰m 嚣p ( x l 铲引 亿 ：三号盟：上 j 纠 鲁 )尸 8 第二章m o n t ec a r l o 方法 图2 2 分布陡峭的被积函数 2 2 3 m o n t ec a r l o 方法在统计物理学中的应用 统计物理学的一个基本任务就是从哈密顿量日出发计算各种宏观物理量。设 有可观测量a ( z ) ，z 是相空间的矢量，是系统处于某状态下一组变量的集合，如 z = ( 盯l ，盯2 ，o n ) 。可观测量a ( z ) 的平均值可由下式定义 1 ， a ( x ) = 专d xe x p 【- h ( x ) k b t a ( x ) ， j( 2 1 3 ) ， 、一。 7 z = d xe x p 【- h ( x ) k b 卅 其中b 是玻尔兹曼常数。系统处在某一状态z 的概率p ( z ) = e x p 一h ( x ) k b t z 。 如果仿照引言中所述的简单抽样方法，随机选取相空间点z f ，然后进行统计平均求 得物理量( 如根据公式( 1 2 ) 计算磁化强度) ，所得结果与理论计算相差很大。原因在于这 样抽样出来的自旋是均匀分布的，其构成的位形也是均匀分布，这样必然得到磁化强度 趋于0 的结果，显然与热力学理论是不符的。在热力学平衡状态，概率函数p ( z ) 分布比 较陡峭，因此对状态z 的抽样应该采用重要性抽样。式( 2 1 3 ) 写成 m e x p - h ( x t ) k b 卅a ( x 1 ) p ( ：r , 1 ) a ( z ) = 曼百一一 ( 2 1 4 ) ， e x p - h ( x 1 ) k b t p ( x t ) 1 = 1 一种简单而且自然的选择是 p ( x 1 ) o ( e x p 【- h ( x 1 ) k b t ， ( 2 1 5 ) 9 第二章m o n t ec a r l o 方法 这时玻尔兹曼因子完全相消，公式( 2 1 4 ) 化为 硐= 丽1 喜m 删 ( 2 1 6 ) 对于已知分布的抽样，可采用的方法有：舍选抽样法、复合抽样法、变换抽样法等。而为 了实现( 2 1 5 ) 式的重要性抽样，直接抽样是很困难的，m e t r o p o l i s 等人提出：不用随机选 取毫不相关的各种状态 z f ，而是构造一个m a r k o v 过程，过程中后一状态觋，由前一状 态z f 通过一个合适的跃迁概率w ( z l _ z f ，) 得到 3 4 。利用m e t r o p o l i s 的想法，在极限 情况m o o 时，通过m a r k o v 过程进行重要性抽样所产生的分布函数尸( z f ) 趋于热平 衡状态下的分布 聃f ) _ 牙1e x p - 掣1 ( 2 1 7 ) 达到这一点的充分条件是加上细致平衡条件 户备( z z ) w ( z l + 宓f ，) = j ( x l ，) w ( 刃f ，+ c ) ， ( 2 1 8 ) 其中墨( 茁f ) 代表系统处于平衡态时的概率分布函数。 显然细致平衡条件( 2 1 8 ) 式并不能唯一确定跃迁概率w ( z z _ z ，) ，其选择仍然有 很多的自由，两种经常使用的选择是 阶= 瓦1 卜a n h ( 器) = 丢煮e x 鬻p - a h k b t ( 2 1 9 ) = = 一一 i 二工7 ， l1 + () 、。 11 t s1 + e x p ( a h k s t ) 或 fi e x p ( - a h k b t ) ，衄 0 ， w ( 规一研) = ( 2 2 0 ) l 三， 日o ， 其中a h = h ( x l ，) 一h ( a z l ) ，7 - s 是一个标度因子，常取1 。可以证明，按照上面两式所产 生的状态序列的概率分布实际上收敛到( 2 1 7 ) 式的平衡概率分布。 2 2 4 单自旋算法和集团算法 ( 1 ) 单自旋算法 在m e t r o p o l i s 构造的m a r k o v 过程中，每次只翻转单个自旋，称为单自旋翻转算法。 翻转概率如果取( 2 1 9 ) 式，称为g l a u b e r 算法( g l a u b e ra l g o r i t h m ) ；如果取( 2 2 0 ) 式，称为 m e t r o p o l i s 算法( m e t r o p o l i sa l g o r i t h m ) 。 1 0 第二章m o n t ec a r l o 方法 采用m e t r o p o l i s 算法进行模拟的过程如下： 1 指定系统的初始状态劫，并设f = 0 。 2 选择一个自旋，假设其翻转，产生新状态z f ，并计算翻转前后的能量差日。 3 如果日 0 ，则按照概率e x p ( 一a h k b t ) 使z f + 1 = z l ，如果日0 ，则直 接令z z + 1 = z z ，。 4 重复第2 步。 当我们选取自旋的次数达到模拟中的自旋总数时，称为一个蒙特卡罗步( m o n t ec a r l o s t e p ，m c s ) 。 g l a u b e r 算法的实现和m e t r o p o l i s 算法基本一致，只是在第3 步所使用的概率不同。 ( 2 ) 集团算法 在临界温度附近，单自旋翻转这种局域算法遇到了严重的临界慢化问题。为了克服 临界慢化，先后提出s w e n d s e n w a n g 【3 5 和w o l f f 3 6 】翻转算法，这两种算法中都是以 自旋集团作为单位来翻转，故称集团算法。 实现s w e n d s e n w a n g 算法的过程如下( 参考图2 3 ) ： 1 任取系统的初始状态劫，并设2 = 0 。 2 遍历晶格所有最近邻座点，如果相邻座点上的自旋值相同，以1 一e 一2 ，蛔t 的概 率将这两个座点用键相连接；如果不相同，则不连接。 3 根据连键的情况构造集团，任何两个座点之间，如果能够找到一条连接的通路， 则属于一个集团，孤立的单个座点也是一个集团。每一个座点必须属于某一个 集团。在确认了每一个集团后，给每个集团按相同几率取自旋的可能值，并完全 忘掉原来的自旋值。 4 完成一个m c s 。重复2 - 3 步。 三三三三三圭誓三二 i 三二；三三三三三 图2 3s w e n d s e n w a n g 算法示例 w o l f f 算法的实现过程如下： 1 随机选择一个座点i ，称为种子座点。 2 从该座点i 出发，如果相邻座点j 上的自旋值与座点i 上的相同，以1 一e - 2 j b t 的概率将这两个座点用键相连，如果相连则将歹座点加入i 所属集团。 第二章m o n t ec a r l o 方法 3 以刚加入集团的新座点出发，遍历新座点的相邻座点，采用上面的规则将更多 的座点加入集团。 4 重复第3 步，直到再不能生成新的键，从而构成大集团，并翻转集团的自旋； 2 3 m o n t ec a r l o 方法的动力学解释 m o n t ec a r l o 方法是通过构造m a r k o v 过程来实现的。在m a r k o v 过程中产生的位 形之间有很强的关联，这种关联对测量值的精度具有很大的影响。这种关联可以通过 m o n t ec a r l o 过程的动力学解释来理解，该解释的根据来自于描述具有随机过程模型 的主方程 3 7 。这种将关联和时间联系起来的解释不仅对理解精度有用，而且是采用 m o n t ec a r l o 方法模拟动力学过程的理论基础。 m o n t ec a r l o 动力学解释是将“时间”t 与演化过程中产生的位形标号联系起来。 假设在时刻t ，位形盯三( 盯l ，0 2 ，a i ，盯) 在m o n t ec a r l o 过程中出现的概率为 p ( o r ，舌) 。概率分布函数p ( o r ，t ) 满足m a r k o v 过程的主方程 t d p ( o r , t ) 一一w ( o r _ o r ) p ( o r ，亡) + w ( o r ，_ o r ) p ( o r , 亡) ，( 2 2 1 ) 7 0- 其中w ( 仃_ o r 7 ) 代表位形由 的概率，称为翻转概率 8 ，3 8 。公式( 2 2 1 ) 右边第一项是使状态o r 概率减少的项，第 二项是第一项的逆过程，使状态盯增加。当系统处于热平衡时，细致平衡条件保证公 式( 2 2 1 ) 中右边两项相互抵消，从而d p ( o r ，t ) d t = 0 。 需要注意的是，一般说来，m o n t ec a r l o 方法产生的多个位形所对应的时间，同真实 物理系统中演化时间是没有关系的。当一o 。时，这个人为随机动力学过程与实际物 理过程都得到相同的热平衡分布。 1 2 j v 研 隔 电 ， 叻 p 吼 至 化变 第三章 二维i s i n g 系统零温动力学 本章我们采用m o n t ec a r l o 方法来研究二维方晶格上i s i n g 系统的零温g l a u b e r 动力 学，包括系统初态对终态的影响，系统的终态位形，以及系统的能量演化等。 3 1 模型和方法 各向同性铁磁i s i n g 模型的哈密顿量可写为 h = 一f 吼 ( 3 1 ) - j “) 在单自旋翻转算法中，单位时间内自旋翻转_ 一的概率为w ( a h ) ，该概率大小 仅取决于翻转前后的能量差 h = 2 h n ， h n ：壹， 3 2 其中h n 是作用在自旋上的局域磁场，z 是晶格配位数，求和对自旋的最近邻自旋 进行。 细致平衡条件( 2 1 8 ) 式可写成 堡! 全堕： a h k