（光学专业论文）准周期十边形覆盖结构上铁磁模型的蒙特卡罗模拟.pdf

摘要 摘要 覆盖理论是对准晶体结构的一个新的描述，与构造彭罗斯拼图需要 两种模块拼砌不同，覆盖结构是由一种结构单元互相覆盖而形成。圉此， 覆盖理论对准晶体的结构及形成过程赋予了新的思想。 本文研究准周期覆盖结构上铁磁模型的临界性质。首先我们从统计 力学的基本原理出发，讨论适合于所研究体系的蒙特卡罗方法，由此确 定了数值模拟所需要的算法及具体步骤。 其次我们研究边界条件问题。由于准周期结构上不能直接应用周期 边界条件，而自由边界条件对计算结果的影响又很大，特别是当系统比 较小时，因此我们用周期覆盖结构做近似，研究准周期覆盖结构上铁磁 模型的物理性质。借鉴产生周期彭罗斯拼图的方法，我们成功地构造出 卜边形覆盖结构上的周期边界条件。 最后我们分别对i s i n g 模型和p o t t s 模型做蒙特卡罗模拟计算。我们采 用w o l f f 算法和单直方图算法得到了临界温度，利用有效尺寸标度获得了 各个临界指数。结果表明基于十边形覆盖结构上的i s i n g 模型与基于正方 品格上的i s i n g 模型属于同一普适类。对于p o t t s 模型也是如此。 关键词：准晶体，十边形覆盖模型，蒙特卡罗模拟，i s i n g 模型，p o t t s 模 型 竺i ! 垦竺! a bs t r a c t t h ec o v e r i n gt h e o r yi san e wd e s c r i p t i o nt oq u a s i e r y s t a l s d i f f e r e n t f r o mp e n r o s et i l i n g sw i t ht w ou n i t s ，t h ec o v e r i n gp a t t e r ni sc o m p os e do fo n e c o n s t r u c t i v eu n i tw h i c hg i v e san e wi n s i g h ti n t os t r u c t u r ea n df o r m a t i o no f q u as i c r y s t a l s i nt h i st h e s i s ，w ei n v e s t i g a t et h ecr i t i c a lpr o p e r t i e so ff e r r o m a g n e t i c m o d e l so nq u a s i p er i o d i cc o v e r i n gs t r u c t ur e w eb e g i nw i t ht h ep r i n c i p l eso f s t a t i s t i c a im e c h a n i csa n dd i s c us sm o n t ec a r l om e t h o dst h a ta r es u i t a b l ef or t h e s y s t e m s s t u d i e d t h e nt h ea l g o r i t h ma n dn u m e r i c a ls i m u l a t i o n s p r o c e d u r ea r ed e t e r m i n e d i nt h en e x tp l a c e ，w e s t u d yt h ep e r i o d i cb o u n d a r yc o n d i t i o n s t h e p er i o d i cb o u n d a r yc o n d i t i o n sc a nn o tb eus e d d i r e c t l y o nq u a s i p e r i o d i c s t r u c t u r e t h ef r e eb o u n d a r yc o n d i t i o n sw i l l a f f e c tt h es i m u l a t i o nr es u l t s ， es p e c i a l l yi ns m a l ls ys t e m s t h e r e f o r ew es t u d yt h ep h y s i c a lp r o p e r t i eso f f e r r o m a g n e t i cm o d e i sb y p e r i o d i ca p p r o x i m a t i o n us i n gt h em e t h o df o r g e n e r a t i n gp e r i o d i cp e n r o s et i l i n g s ，w ec o n s t r u c tt h ep e r i o d i cb o u n d a r y c o n d i t i o nso nd e c a g o n a lc o v e r i n gs t r u c t u r e f i n a l l y ，w ed o m o n t ec a r l os i m u l a t i o n sf o ris i n gm o d e la n dp o t t s m o d e lo nt h ed e c a g o n a lc o v e r i n gs t r u c t u r e w eg e tt h ecr i t i c a lt e m p e r a t u r e t h r o u g hw o l f fa l g o r i t h ma n ds i n g l eh i s t o g r a mm e t h o d b yu s eo ft h ef i n i t e s i z es c a l i n gw eg e tt h ec r i t i c a le x p o n e n t s i tisf o u n dt h a tt h el s i n gm o d e lo l l t h ed e c a g o n a lc o v e r i n gs t r u c 。t u r eb e l o n g st ot h es a m eu n i v er s a lc l a s sa st h a t o ft h es q u a r el a t t i c e s t h es a m ec o n c l us i o ni so b t a i n e df o rt h ep o t t sm o d e l k e y w o r d s ：q u a s i c r ys t a l ，d e c a g o n a l c o v e r i n gm o d e l ，m o n t e c a r l o s i m u l a t i o n ，i s i n gm o d e l ，p o t t sm o d e l l 华南理工大学 学位论文原创性声明 本人郑重声明：所呈交的论文是本人在导师的指导下独立进 行研究所取得的研究成果。除了文中特别加以标注引用的内容 外，本论文不包含任何其他个人或集体已经发表或撰写的成果作 品。对本文的研究做出重要贡献的个人和集体，均已在文中以明 确方式标明。本人完全意识到本声明的法律后果由本人承担。 作者签名：写伛跌 日期：a f 年多月f t 日 学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规 定，同意学校保留并向国家有关部门或机构送交论文的复印件和 电子版，允许论文被查阅和借阅。本人授权华南理工大学可以将 本学位论文的全部或部分内容编入有关数据库进行检索，可以采 用影印、缩印或扫描等复制手段保存和汇编本学位论文。 保密口，在年解密后适用本授权书。 本学位论文属于 不保密q 。 ( 请在以上相应方框内打“”) 作者签名： 导师签名： 乃位玖 懈 日期：枷5 年g 月陟日 日期：埘年舌月, w h 第一章绪论 1 。1 准晶体概述 第一章绪论 在传统的固体物理学中，固态物质可分为晶体和非晶体两大类。晶 体可以产生明锐的布拉格衍射斑纹，这是由于原子在晶体中的周期排列 使得衍射束具有高度的相干性，导致出现很强的衍射峰。而非晶体由于 原子排列是杂乱的，因此只能产生弥散的衍射谱。晶体的周期性使得在 晶体中不能存在五次及大于六次的旋转对称性。 但是在1 9 8 4 年，美国国家标准局s h e c h t m a n 等人在研究合金时， 意外地发现了一种a 1 m n 合金具有五次对称性的衍射图，从而揭开了准 晶体研究的序幕。 准晶体是一种具有长程准周期平移序和长程取向序物质结构的晶 体，它的研究首先从实验开始，对各类准晶体材料的制备测试和结构分 析已比较成熟，对准晶体的力学、热学、光学、电学等物理性质的实验 研究也取得了很大的进展。而理论研究方面则相对落后，其主要原因在 于：长程平移序的缺失使得解析解难以得到，而简单的构形平均又不能 反映准周期的特征，因而数值模拟计算便成为主要的研究手段，但数值 模拟有计算机容量及速度上的限制，所以在大多数的情况下只能得出一 些近似的结果。 准晶体的结构模型主要有两种，其中彭罗斯拼图( p e n r o s et i l i n gs ) ( 2 】是由两个不同的单元按照一定的拼砌规则，在平面拼砌而成的二维准 周期结构。彭罗斯拼图早在准晶体发现前已经由英国数学家彭罗斯发现 了，它的傅立叶变换具有五次对称性。鉴于彭罗斯拼图和准晶体的相似 性。它很快就被作为准晶体的结构模型。彭罗斯模型成功地解释了一些 准晶体的结构和物理性质。然而，对这个理论的雁确性也一直存在疑问： 构造彭罗斯拼图至少需要两块拼块，为什么一种物质中要有多种藿复单 元? 彭罗斯拼图的拼砌规则在准晶体的形成过程中是如何实现的? 带着疑问，德国青年数学家g u m m e l t 对彭罗斯拼图做了进一步的研 究，提出了十边形覆盖模型( d e c a g o n a lc o v e r i n gm o d e l ) 【3 】，它不同于 彭罗斯拼图，只由一个单元( 双色十边形) 按_ 定的规则覆盖就能形成：维 五次准周期结构，且同样具有五次对称性。尽管它跟彭罗斯拼图存在着 定的对应关系，但是它具有不同的物理内涵：在十边形覆盖模型中只 华南理工大学硕十学位论文 存在着一种重复单元，可称之为准原胞：相邻单元的覆盖是近程操作， 只要求双色十边形的同色部分重叠；重叠的部分已经在原有的十边形单 元中存在，可以作为新生长出的f 边形的单元核。一1 些实验和理论研究 对十边形覆盖模型给以支持t 4 。6 l 。 彭罗斯拼图和十边形覆盖模型都是二维五次对称性的。从对称性束 看，还有八次、卜次、十二次等高次对称性的准晶1 7 。从维数来看，还 有一维准晶和三维准晶。一维准晶最简单，研究已经很成熟了；而三维 准晶由于结构过于复杂，且由于计算机的速度容量的限制，理论研究较 少。因此目前研究得比较多的是二维准晶。 对于准晶体物理性质，虽然无法解析求解，但可以用很多近似的方 法处理，这是很多科学家常用的方法。用周期彭罗斯拼图( p e r i o d i c p e nr o s et i l i n g s ) f 8 】近似代替彭罗斯拼图，计算准晶体的物理性质就是这 样的一种思想，而且已经取得了一定的成果。在本文中我们将此方法应 用予十边形覆盖模型上，构造出周期单元尺度不同的周期十边形覆盖结 构( p e r i o d i cd e c a g o n a lc o v e r i n gs t r u c t u r e ) ，进而研究其物理性质。 1 2 铁磁模型的临界性质 相变和临界现象是凝聚态物理中一个很重要的课题【9 12 1 ，这是因为 临界区域是物质不同相的过渡区域，在这一区域的物理性质丰富且有深 刻的意义。在铁磁模型中，研究得最多的是i s i n g 模型【1 l ，l3 了。丁f 方晶 格i s i n g 模型是可解析求解的，已精确求得其临界温度 t = 2 1 1 0 9 ( 1 + 、2 ) = 2 2 6 9 ，而其i 临界指数有：u = 1 ，y = 1 号，叩= 告，声= 【i4 1 。 自从准晶体被发现以来，得到了众多学者的关注，但是由于其准周 期性，使得无法求得建立于准周期结构上的is i n g 模型的解析解，只能通 过数值模拟的方法来近似计算其物理性质。在周期边界条件f ，已计算 得到彭罗斯拼图is i n g 模型的临界温度= 2 4 0 1 + 0 0 0 5 ，其临界指数”，r 在误差允许范围内跟正方晶格i s i n g 模型的临界是相同的，这说明两者属 于同一普适类f 15 】。而已知建立于非周期结构上的铁磁模型可能有彳；同的 普适类，因此有必要研究基于覆盖结构上的铁磁模型的l 临界性质。在f - 边形覆盖结构上的i s i n g 模型的临界性质是本文的主要研究对象。p o r t s 模型1 1 1 ，13 ，16 l 是另一个得到广泛研究的铁磁模型。我们也对基于十边形覆 盖结构的临界性质做详细的研究。 第一章绪论 1 3 本论文的主要内容 本论文研究建立于周期十边形覆盖结构上i s i n g 模型和p o t t s 模型的 临界性质，内容主要包括在以下各章中： 第二章介绍。些背景知识：包括统计力学原理，要研究的物理模型 一一l s i n g 模型和p o t t s 模型，磁体的临界性质以及用于计算临界指数的 有限尺寸标度。 第三章阐述在模拟计算中要使用的方法一一蒙特卡罗方法。包括蒙 特卡罗方法的原理、以m e t r o p o l is 算法为例说明的模拟过程、专门针对 l 临界点模拟的w o l f f 算法、针对蒙特卡罗方法的数据处理方法一一单直方 图方法。 第四章研究如何把周期性边界条件推广到十边形覆盖模型中去。这 将从彭罗斯拼图的构造出发，引出广义对偶法，利用广义对偶法产生周 期彭罗斯拼图，在这样的拼图中两种菱形并不完全满足彭罗斯拼砌规则： 接着研究彭罗斯拼图和十边形覆盖模型的具体对应关系；最后从周期彭 罗斯拼图得到周期十边形覆盖结构。 第五章是在前面三章建立的原理和方法的基础上对铁磁模型做蒙特 卡罗模拟。首先选取模拟的系统，它们由数目及尺寸不同的周期十边形 覆盖模型中的周期单元组成；然后确定要模拟系统的平衡时间和相关时 间。为了提高模拟的效率，采用了针对临界区域的w o l f f 算法和单直方图 方法，得到系统的临界温度和临界指数。从临界指数的模拟结果，可以 判断基于十边形覆盖结构的i s i n g 模型和基于正方晶格的i s i n g 模型是属 于同一普适类的；其临界指数是符合标度假设的。对于p o t t s 模型也是如 此。 第二章基础理论 第二章基础理论 本论文的工作是应用蒙特卡罗方法【l3 】来模拟计算准晶体的物理性 质，而蒙特卡罗方法是以为统计力学为基础的，因此我们首先需要了解 统计力学的一般原理。 2 1 统计力学背景 期望值 统计力学主要用统计的方法研究凝聚态物质系统的物理性质。凝聚 态物质系统是一个多体系统，包括有巨大数目的单体。每个单体都有各 自的动力学方程。因而对凝聚态物质系统来说，其动力学方程的数目是 惊人的，这使得无法严格地做数学求解。表面看来，似乎很难知道一个 凝聚态物质系统的性质了，但现实情况是，系统的行为往往是可预测的， 比如气体的压强和温度之间的关系服从简单的定律。统计力学绕过了求 解多体系统的动力学方程组的难题，它认为系统的性质是对系统各个微 观状态做统计平均的结果，可以用统计的方式来计算大系统的宏观性质。 比如我们平常测量某个凝聚态系统，因为系统微观状态的变化非常快， 在测量过程中，系统已经历过了很多个状态，测量的结果是这些状态的 统计平均结果。 最典型的是哈密顿系统：系统有一定的哈密顿量，处于某些微观状 态中( 或者说系统处于具有一定哈密顿量的任何状态中) 。除了哈密顿系 统还有热源，热源和哈密顿系统之间相互交换能量。热源能推动系统从 一能量状态到另一能量状态。 假设系统处于一定微观状态1 中，p ( - 秒) 是从状态到”的迁移率， 一股认为迁移率与时间无关；w ，( f ) 表示在t 时刻系统处于状态k t 的几率， 那么我们可以得到一个主方程： d w u ( t ) ： ( f ) p ( v - - 4 t ) 一( f ) p ( , u - 9v )】。dt 1 ”、 。- ” ”。 当然，几率w u ( t ) 必须满足归一化条件： ( f ) = - 。 ( 2 1 ) 假设要测量物理量q ，可以通过求其期望值来得到。在知道物理量q 在各 状态2 下的几率w u ( t ) ，以及其取值线，那么可以得到其期望值为： = z q u w u ( t ) 。( 2 3 ) n 华南理工大学硕士学位论文 这个期望值尽管不一定跟实际观察值完全符合，但其包含了物理量q 的 重要信息，在一定程度上可代表物理量q 的实际观察值，一般实际观察 值就在期望值附近。 物理量q 的期望值 可以理解为系综平均的结果：有许多相同而 独立的系统，各自拥有自己的热源，各个系统在自己热源的作用下，独 立而快速地从一状态变到另一状态，对各个系统的观察值q 做统计平均 就得到其期望值 。只要取的系统足够多， 会是一个很好的估计 值。但是在实际中我们无法拥有许多模一样而独立的系统。特别是在 数值计算中，由于内存空间的有限，这种方法是行不通的。 另一种不够严格的理解方法是：对于一个系统，每隔一定时间测量 并得到其观察值q ，然后再对所有观察值q 取平均而得到其期望值 。 其不严格之处在于两个时间间隔短的状态存在一定的相关性，而我们整 个测量过程的时间是有限的，在有限长的时间内可能无法抽取到足够多 的有代表性的状态。 平衡态 我们来看看方程( 2 一1 ) ，当方程右边为零时表示系统处于a 状态的 几率不变了。如果各个状态的几率都不变了，那就认为系统平衡了。一 般认为由方程( 2 1 ) 和( 2 2 ) 决定的系统，最终都会达到平衡状态。 达到平衡状态的w 。( f ) 称为平衡占有率，表示为： 以2 烛( f j 。 ( 2 4 ) 吉稚斯于1 9 0 2 年证明了对于与热源处于平衡状态的系统，假设温度为了1 ， 平衡占有率满足玻耳兹曼分布： 1目 儿= - b e 7 ”。 ( 2 5 ) 厶 这里e 。就是在状态时系统的总能量；k 为玻耳兹曼常数；z 是配分函数， 表示为： z = p 一= p 一肫。 ( 2 6 ) 一一 。、7 f 其中= ( 灯) ，称为反温度。如果能够知道配分函数随温度和影响系统的 其它参数的变化，那么就能够知道系统的所有宏观行为。 从式( 2 3 ) ，( 2 4 ) 和( 2 5 ) 可得，对于一个处于平衡态的系 统其物理量q 的期望值为： = g 砟= i 1 q p 邓丘。 ( 2 7 ) “厶“ 例如，内能u 的期望值为： 6 第二章基础理论 u = 专莓q 。饵= i 1 两b z a 矿l n z 。( 2 - - 8 ) 比热c 是系统内能的微分为： c = 券硼2 等邶2 等。 cz 吲 同时，在热力学中，比热c 和熵s 的关系为： c = 丁万0 5 s = 一喏。 c ： 从( 2 9 ) 、( 2 1 0 ) 两个式子，可得出的函数，对积分，可得到 熵s 为 沣t 卢等z 。 从( 2 8 ) 、( 2 1 1 ) 两式子可得到系统的自由能f 的表达式 f = u 一碍= - k t l l l z 。 从( 2 8 ) 、( 2 1 1 ) 两式子可得系统内能u 和自由能_ f 的关系 r 2 警一r 2 警。 ( 2 12 ) ( 2 一l3 ) 通过上面的式子，可以从配分函数直接计算内能，自由能，比热和 熵了，除此之外，还可以通过自由能计算其它影响系统的参数。在经典 热力学中，和系统相互作用的外场或者外部约束都有其共轭变量，以表 征该外部作用对系统扰动的响应。比如，一个理想气体系统的体积y 变 化导致了其压强p 的变化，理想气体的体积v 和压强p 是一对共轭变量。 磁体的磁化强度m 的变化是对作用在磁体上的磁场曰的响应，磁体的磁 化强度肼和作用于其上的磁场强度曰是一对共轭变量。自由能对某个外 部作用的微分可得其共轭变量。如： a f p 2 一万 m = 誓。 ( 2 14 ) ( 2 15 ) 涨落，相关性 统计力学除了能得到系统的宏观性质外，还能得到经典热力学所无 法得到的一些性质，比如可观测量的涨落。一般来说，我们除了计算可 观测量的平均值外，还有必要计算其标准差，因为它能衡量平均值的近 似程度。以系统内能为例，其标准差为； ( ( e 一 ) 2 ) = 一 2 。 ( 2 16 ) 而 7 华南理工大学硕士学位论文 = 三zy ，f 2 一一辑= 上丝a f t 2 z 。 所以 ( ( 趴酬= i 1 矿0 2 z c 舄2 = 等。 用式( 2 9 ) 消去二次微分，可得： ( 趾纠2 # 斋。 式( 2 1 9 ) 有几个意义：第一，它把内能的涨落以比热的形式给 出，换句话说就是能用经典热力学的物理量表示涨落，虽然涨落无法在 经典热力学中求得，因为涨落依赖于微观细节，而经典热力学不研究微 观细节。第二，对于一个大系统来说，可以看出涨落是极其微小的。第 三，比热是广延量，与系统体积y 成正比，而内能的涨落实际是式( 2 一 1 9 ) 的开方，也就是说与y 成正比；而我们知道内能本身与体积成f 比， 这就是说内能涨落的相对大小与咖成反比，即系统越大，其内能的相对 涨落越小。这一结论可以推广到其它可观测量上。可以想象当系统足够 小的时候，涨落就不能忽略了。系统小到不能忽略涨落的极限为热力学 极限。 除了内能的涨落，还可以计算其它热力学量的涨落。_ | 二面已经提到， 对应于系统的外场或者外部约束都有一个共轭变量。假设外场为f ，其对 应的共轭变量为x 。，在系统的内能e 包含有多个一x 。l ，n 表示各种外场 或者外部约束及其共轭变量，其中的各个e 是相互独立的。因此x 。町由 下面公式给出，这里为了简单，把它写为x ，其对应的外部变量写为y ， 则： = 去莓x e 叩巳2 瓦1 面0 p l - 7 e 一盹 up lo l 也可以写成： ， 1a l n z b f 口a ya y 。 那么w 观测量x 的涨落： 1 0 3 4 k 时，铁是顺磁 质，只有在加上磁场日的情况下，铁才被磁化，且磁化强度跟磁场h 成 下比。而当温度下降到1 0 3 4 k 以下时，则即使无外加磁场的情况f ，铁仍 带有磁性，且其磁化强度跟加于其上的磁场不是线性关系。铁磁体的相 变过程中单自旋磁化强度m 是连续的，而其一阶导数则是不连续的。i s i n g 模型作为一个磁体模型，其相变跟铁的相变是同一性质的。 相变过程中还有其它的特点，比如在临界温度附近l s i n g 系统的波动 很剧烈，而且越靠近临界温度，这种波动越剧烈。这将导致临界点附近， 系统一些性质以幂率形式变化。对于某一类物质系统，这些幂率中的指 数是常数，称为临界指数。通过临界指数可以知道系统的重要属性，比 如其相交的级数。 本论文研究i s i n g 模型和p o t t s 模型在准晶格上的临界性质，因此在此 只介绍磁体的i 晦界指数。在后面对讨论中，没有特别说明，也都是以这 两个磁体模型为例的。 先定义约化温度： r = 等， cz 吲， 其中z 表示临界温度，用约化温度来表达系统在临界点的性质比较方便。 华南理工大学硕士学位论文 磁体的性质可以表示成温度和外加磁场的h 的函数， 磁化强度m 和磁化率z 的函数关系如下： m ( o ，t ) ( 一f ) 4 当t _ 旷时 m ( h ，r ) h y 8 当f - - - ) o h 豕o ，t ) t 1 当t _ 0 + 时 z 一 。 ( 2 3 4 ) 司以用关联函数d 2 ( i ，) 来衡量两个不同格点间的相关程度，相关程度随 着两点之间的距离增加而减小，因此也可以用相关长度毒( 0 ，t ) 来衡量相关 性，它表示相关函数降到一个比较小的值时两点的距离。相关长度在临 界点的行为是： 善( o ，丁) t “当t _ 0 + 时 ( 2 3 5 a ) 善( o ，r ) ( 一r ) “ 当t _ o 一时( 2 3 5 b ) 其中”，秒为临界指数。把关联函数q 2 ( f ，j ) 在临界点的相关性表示为g ( r ) ， r 为两点之间的距离。关联函数的临界行为表示为： 占( r ) r “2 一， ( 2 3 6 ) 其中d 为系统的维数，玎为临界指数。 标度假设与临界指数普适性 自然地，人们希望寻找到简单表达式来描述上面的临界点附近的热 力学函数。w i d o w 、d o m b 和h u n t e r 认为某些热力学函数可能是相似的。 g r i f f i t hs 特别指出磁场强度h 和m “，与t 的函数关系是相似的。既然磁场 强度h 是磁化强度 f 的奇函数，则在临界温度r 附近有： h k t 。= m i m i “1 h , l m i 一1 7 9 ) ， ( 2 3 7 ) 其中的和j 是待定的常数，k 是玻耳兹曼常数，而矗。( 工) 是无量纲的标度 函数。如图2 - 2 是正方晶格i s i n g 模型的标度函数。整个函数是连续的， 在一k x 7 1 l i r a g ( t , l 卜 1 ：r ， ( 2 4 4 ) 而在有限系统中，其可表示为： g ( t ，l ) = g ( t l r ) ， ( 2 4 5 ) 其中t 为约化温度。在临界温度r 下，t = 0 ，这时b i n d er 四阶累积量跟系统 大小无关，通过对不同尺寸的系统的模拟，求得相应的b i n d er 四阶累积量 相对于温度的关系图，其交点就是临界温度所在了。 当然模拟必须建立在相同几何结构上的，因为尽管对于不同的几何 结构临界指数可能是一样的，但是系统的性质是跟系统的几何结构有关 的。临界指数只能反映系统性质的一部分。 再进一步，对b i n d e r 四阶累积量求导，可得到临界点处的导数为： 把= 等 通过上式可以求得临界指数”。 其它的系统性质也有相应的有限尺寸标度。 m ( t ，l ) m 2 ( i ，) ”2 z ( i ，l ) ” c ( l ，l ) 一 ( 2 4 6 ) ( 2 4 7 a ) ( 2 4 7 b ) ( 2 4 7 c ) ( 2 4 7 d ) 如果对上面的式子左右两边取对数，则可以用直线拟合的方法求得式中 的临界指数。 第三章蒙特卡罗方法 第三章蒙特卡罗方法 从2 1 节中知道，计算一个凝聚态系统的可观测量的期望值可以通过 求其配分函数的方法来得到，但是这种方法在数值计算中是行不通的。 比如对于一个格点数为的i s i n g 系统，其状态数是2 ”，对于一个比较 大的系统，其状态总数是惊人的。一般我们只能抽取有限的状态做统计 平均，这种方法是不严格的，因为抽取到的状态可能没有代表性。但是 如果让模拟过程中抽到某一状态的几率反映其理论上存在的几率，结果 就不同了，蒙特卡罗方法为了实现这个目的的一个方法。 本章将讲述对系统的热平衡态的蒙特卡罗模拟。首先是其原理，其 次是具体的算法，最后是具体的数据处理方法。 3 1 平衡态蒙特卡罗方法原理 估算子 对热平衡系统做模拟计算可以求可观测量的期望值，比如内能、比 热、磁化强度、磁化率等。比如求可观测量q 的期望值，理想的方法是对 系统的所有状态以玻耳兹曼分布率为权重来求平均： q x 挈妥，( 3 - - 1 ) 但这种方法只有在系统极小时才有可能。对于大的系统，只能抽取- 部 分的状态来求平均。抽取状态是以某一几率分布见为基础的，这个跟模 拟的方法有关了。假如抽取了状态子集f 麒心 ，则： ：喾，( 3 - - 2 ) 线称为q 的估算子。 一般我们选择的抽取几率是对各个状态都是相等的，所以式( 3 2 ) 可简化为： 线：圣趋： 很明显上式的结果是很不准, ”, 确e - e 的e j ，因为我们只能抽取一小部分的( 系3 统- - 3 状) 态，而且可能抽取到一些权重很小的状态。但是，如果我们能预先知道 某些状态在式( 3 一1 ) 中起重要作用，而且让在抽取过程中抽中这些状 态的几率比较大，那么计算的结果就很准确了。 华南理1 ：大学硕十学位论文 重要性抽样 重要性抽样就是把那些具有重要贡献的状态抽取出来的抽样方法。 我们知道，系统的热平衡态占有率是玻耳兹曼分布率，如果以玻耳兹曼 分布率作为抽取几率来抽取状态，这是最好不过的，这样式( 3 2 ) 就 司简化为： o e u e i ( 0 。 以式( 3 1 2 ) 为选择率，式( 3 一1 4 ) 为接受率的计算就是著名的m e t r o p o l i s 算法了。 平衡时间 我们要计算处于平衡态的系统的性质，但模拟总是从一个初始状态 开始，而初始状态并不就是平衡态。从初始状态到平衡状态需要经历 段时间，这段时间就称为平衡时间，以瓦。表示。 那么如何看出平梃时间。一般每走一个蒙特卡罗步( 就是系统的每 个格点都有一次反转的机会) 就计算一次系统某可观测量的值，得到一 系列值后就可以画一个该可观测量与蒙特卡罗步之间的关系图了，如图 3 1 ，从图中就可以估计出系统的平衡时间了。如果看出该可观测量经过 一定的时间后，测量值比较稳定了，就说明系统可能处于平衡状态了。 但不能绝对肯定，因为如果系统处于亚稳态区域也有可能出现这样的情 况。为了减小误判的几率，一般可以从多个初始状态开始，画多个关系 图，如果各个关系图的最终趋势相同，那基本上就可看出平衡时间了。 对于无外场的i s i n g 系统，系统的平衡时间跟两个因素有关。一是温 度，越靠近临界温度，系统的平衡时间越长。另一方面是系统的大小， 系统越大平衡时问越长。 图3 1平衡时间 2 0 量瞢每cq! 第三章蒙特卡罗方法 相关时间 当系统达到平衡后，还必须注意一点，那就是相关时间f 。简单的说， 相关时间就是系统从一个状态到另一个与原状态明显不同的状态”所 需要经历的时间。 我们知道，新状态是在旧状态的基础上产生的，两者是有一定的关 联的，在时间上相近的状态有一定的相关性。一般来说，时间上相距越 远，相关性越弱。蒙特卡罗模拟过程中，只抽取一极小部分的状态出来 计算，比如对于一个有1 0 0 0 个格点的系统，其系统总状态数为2 1 ”，这 足一个很大的数目，而一般我们抽取的状态数只有十万个，百万个，相 比较而言是很小的，所以抽取出来的状态必须非常有代表性，而不是一 系列相关的状态，不然计算的结果就不可靠了。一般是每隔一段足够长 的时间抽取一个状态。 相关时间f 可以通过自相关函数来衡量，自相关函数是这样定义的： r ( t ) = i 由，【m ( f ) 一 【m ( f + f ) 一 】， ( 3 15 ) 其中m ( t ) 为单自旋磁化强度对时间的连续函数。若m ( t ) 是离散的，则： 心，2 去和撇卜去和勺击踟， ( 3 一1 6 ) 其中f 是两次测量的时间间隔，f 一是从测量开始到结束的整个时间。一般 自相关函数正比予p 一以，相关时间f 就是自相关函数降至z 的时问。相隔 一个相关时间的两状态仍有一定的相关性，为了使相关性更小些，可以 取更长的抽样时间间隔，一般两倍相关时间足够了。为了计算相关时间 的方便，一般是采用积分的形式： f e - t l 出= 熟 这样比直接用r ( f ) 去拟合指数函数求f 更好些。 m e t r o p o l i s 计算的相关时间是跟温度相关的，一般越靠近临界点，相 关时间越长，这是m e t r o p o l i s 算法的缺点，针对这个缺点，发展了w o l f f 算法。 m e t r o p o i is 算法的具体步骤 至此我们已经介绍完了蒙特卡罗模拟过程中所需要注意的问题了， 现在来看看模拟过程的具体步骤。一个蒙特卡罗步包括如下过程： 在晶格中随机选取一点； 计算将该格点反转后系统能量的变换胡； 华南理工大学硕士学位论文 按照m e t r o p o l i s 算法的接受率判断是否接受反转： 假如系统有个格点，重复步骤次。 一个蒙特卡罗步就是要让所有格点都有一个反转的机会。如果随机数的 质量够好，做次随机选取格点跟直接对格点进行逐一选取基本上是一 样的，用逐一选取能节省产生随机数的计算时间。 以上是一个蒙特卡罗步的具体内容，把它当作一个功能单位来描述 整个m e t r o p o l is 算法： i 对系统进行初始化，一般是让晶格能量最低，即所有格点的自旋 具有同一取值； i i 做次蒙特卡罗步，以让系统达到平衡状态( 乙就是前面所讲的 平衡时间，必须提前算好) ； i i i 每隔若干个蒙特卡罗步抽取一个状态，至于要隔多少个蒙特卡罗 步要看相关时间的多少了： i v 用抽取的状态进行计算，如何计算在后面会讲到： v 重复步骤i 、i i 、i i i 、i v 若干次，从几万到上百万不等，越多计 算结果越可靠。 这里是以m e t r o p o l i s 算法为例讲的，对于后面的w o l f f 算法、热浴算法， 基本步骤是相同的。 3 3w o l f f 算法 临界点的相关时间 m e t r o p o l i s 算法在模拟远离i 临界状态的系统是很准确的，但在临界温 度附近有很大的误差。这里面有两个原因，一是临界波动增大了，一是 临界温度附近相关时间变长了。临界波动增大是i s i n g 模型的固有性质， 无法改变，而相关时间则可通过改良算法来改善。下面先介绍系统临界 点相关时间的特点。 处于平衡态的i s i n g 模型会形成一簇簇具有相同指向的格点集，格点 集的典型尺寸称为相关长度，用手表示。相关长度孝随着温度趋近于z 而 趋近于无穷。跟第二章中一样以约化温度t 为变量： 扭等。 相关长度在临界点附近的发散程度可以表示为 孝圹， 其中p 是一个正常数，称为临界指数。 第三章蒙特忙罗方法 相关时间也有类似的发散效应，但相关时问跟算法有关，其发散形 式为： f 一旷7 ， ( 3 2 0 ) 其中z 称为动力学指数，对于不同算法，可以有不同的取值。上式显示出 了越靠近临界温度，相关时间将越长，直至无穷，这使得在临界温度的 模拟运算很慢。 但实际上相关时间并不会趋于无穷，因为我们采用的是有限系统。 把式( 3 1 9 ) 和( 3 2 0 ) 合并起来，得： 彳一f 。 ( 3 2 1 ) 而对于有限系统，假设其边长为l ( 对于二维准晶体，可以用晶格点数开 方得到) ，则孝的最大值只能为l ，所以在临界点附近的相关时间可以表 示为： f 口。( 3 2 2 ) 研究结果表明对于m e t r o p o l i s 算法，z 约为2 ，对于一个较大的系统，上 式的相关时间是无法忍受的。为了模拟临界状态，必须改变算法减小动 力学指数z 。 w o l f f 算法 由于在相变点附近，相关长度很大，以i s i n g 模型为例，系统中存在 很大的自旋方向相同的区域。在m e t r o p o l i s 算法中，我们使用的是单自旋 反转，反转这些区域时显得很慢，这是导致m e t r o p o l i s 算法动力学指数大 的原因。为了加快系统变化的速度，我们采用簇反转的算法，就是把具 有相同自旋方向的格点所组成的簇做整体反转。w o l f f 算法 2 0 1 就是这样 的一个算法。 在晶格中随机选取一个格点，然后看其邻近格点是否与其自旋方向 相同，如果是，再看这些邻近格点的邻近格点如何，如此重复，直至找 出整个簇。当然，不能把整个簇反转，反转的数目要依赖于温度。因为 随着温度的下降，簇的尺寸变得越来越大，所以要让反转的数目随着温 度的下降而递增。在w o l f f 算法中是这样实现的：不让所有与第一个选取 格点自旋相同的格点加入到簇中，而是以一定的几率见。加入，而以。随 着温度下降丽上升，最后把这样形成的簇以依赖于反转所消耗的能量的 接受率做反转。 如图3 2 中的两系统状态肛和口，因为t 和 两状态中所选取的簇是 一样的，若以上面所描述的方法来选取簇，则两簇被选中的几率是一样 的。不同的是所要挣破的束缚不一样。要挣破的束缚就是那些与簇同向 华南理t 大学硕十学位论文 但没被加入到簇中的临近自旋。有一个同向而不被加入到簇中的几率是 1 n 。，设在中有n 个这样的自旋，在口中有m 个这样的自旋，则选择率 g ( 。 ) 为( 1 - p o a d ) ”，选择率g ( v _ ) 为( 1 一儿。) “，则细致平衡条件为： 墨! 壁三塑垒! 岂2 型： g ( v - - f 1 ) a ( u 斗)( 1 一p 。“) j a i ( f l 而- - - ) v ) = e 一4 ( 毛q ) ， 而在外加磁场为零的情况下，系统能量的变化为 r 一瓯= 2 j (