式给出的平衡概率P。(x2)2只要N。/N,小于平衡概率之比, 我们就有ΔN→。>0,即比值N。/N,向着平衡概率的比值增 大;反之,如果N。/Nr大于“平衡比值”,那么△N→,<0, 从而N/N,再度减小到正确的平衡比值。于是当l→∞时将渐 近地到达一个定常态分布,它的N。/N,之值正是平衡分布 (2.1.36)式所要求的。代替考虑同时有多个 Marko链,我们 也可以等价地把一个很长的 Markov链切成等长的若干段,而对 链的相继各段应用同一论据 现在我们来讨论这个问题:变迁x4→xt,在实际做时应当怎 样取?在原则上,这个变迁的选法有很大的自由,唯一的限制是 “先验概率”W。*=0(x1→x…)应当对称这个条件,即 WF=0(xt→x1)=W 并且由此得出的在有能量变化8时的跃迁概率τ,W(x1→x1,) 之值,应当以足够的频度显著异于0和1。因此最常取的变迁 是,只在一个(或几个)自由度上发生变化因为如果我们同时 改变N》1个自由度,那么(2.1.39)式中的8/BT就将是 N(E/kT)的量级,其中E是能量标尺(例如,对于磁哈密顿量 (2.1.1)-(1.3),E=j),由于感兴趣的温度是使E/kBT的数 量级为1的温度,那么对于几乎每一个涉及N>1个自由度的变 迁,如果它要消耗能量,它的跃迁概率都是极小的,因而绝大部 分变迁的企图都根本不会执行;系统会“粘死”在它前一个位形 上。这在大多数情况下显然导致一个不实用的算法。克服这一困 难的一个办法是联合的蒙特卡罗-郎之万算法t2 2.1.7再谈模型和算法 图2.4示出对统计力学中研究的若干模型常用的一些位形变 迁。对于 Ising模型,最常用的算法是单个自旋反转算法和自旋 交换算法(图2.4(a)(b))。注意单个自旋反转算法显然不会
使系统的总磁化强度为一不变量,而自旋交换算法则使总磁化强 度不变。因此,这些算法对应于实现不同的热力学系综:图24(a) 实现一个“巨正则”系综,这时温度T和磁场强度H是独立给 定的热力学量,需要计算的是其共轭热力学量(例如,磁化强度 <Mr是和磁场强度H共轭的量),而图2.4(b)则实现一个 “正则系综”,这时温度T和磁化强度M是独立给定的热力学量 (而磁场强度<H>r是我们想要计算的共轭因变量) 我们把 Ising模型的(T,H)系综叫做“巨正则”系综,而 把(T,M)系综叫做“正则”系综,这用的是 Ising模型的格气模 型解释的语言,在这种解释中,自旋变量S;通过p;=(1-S;)/2 被重新解释为局部密度p(=0,1)。于是<M>r同平均密度<0;>r 相联系:<Mr=1-2<P;;而H则同可能占据点格座点的粒子 的化学势相联系。 已知在热力学极限N→∞下,统计力学中不同的系综给出等 效的结果。于是对统计系综从而对与之相联系的算法的选择似乎 只是一个选哪一种更方便的问题。但是,各种系综中的有限尺寸 效应很不相同,而且在一次模拟中趋于平衡位形的“速率”也可能 很不相同,因此选择适当的统计系综有时是一个复杂的问题。当 我们用“速率”这个词时,我们心中已想到了蒙特卡罗过程的动力 学解释20:于是图2.4(a)实现 Glauber2110动态 Ising模 型,它是一个单纯弛豫模型,没有任何守恒定律;而图2.4(b)则 实现川崎恭二212的动态 Ising模型,其磁化强度守恒。 对于具有连续自由度的模型,例如ⅹY模型(21.2)或 Heisenberg磁体(2.1.3),还有对于流体的模型(图2.4(c)(d)) 不是完全随机地、而是在其先前的值周围的一个区间内选取一个 粒子的新自由度(例如图2.4(c)中的角度y,或图2.4(d) 中的位置x,v;),常常是合适的。然后对这个区间的大小加 以调整,使对图24中考虑的试验变动的平均接受率不至于太小
t→ 4t11=111t 2△ =g+△(25-1) 2△y{2→; (d) 2△xx,"x+△x(2-1),yy,+△y(2-1) 八 八 (e) 中一中 utr 图2.4统计力学的一些标准模型的蒙特卡罗模拟中常用的位形变迁 xx1的例子。(a)单个自旋反转 Ising模型(按照动力学解 模型(按照动力学解释,这是川崎动态 Ising模型),()的 释,这是 Glauber动态 Ising模型);(b)最近邻交换Isin 个在[0,1]区间均匀分布的随机数,XY模型的算法的两种变型 左图,描写自旋的新方向的角度卩;完全随机地挑选,右图,φ;从 原来的方向P两侧的区间[;-△甲,卯;+△]中选出,(d)二维 流体中一个原子的坐标变迁,从旧位置(x;,3)变到一个新位置, 这个新位置均匀分布在旧位置周围的一个大小为(2△x)(2△y)的 方形中;(e)在一个单座点位势V(φ)中运动的一个粒子的变 迁,从旧位置ψ变到新位置中;(f在聚合物的点格模型和非 点格模型中使用的变迁(被移动的键用波纹线表示):左图,方形 点格上的蠕行算法;中图,与聚合物动力学的 Rouse模型有关的 个动力学算法;右图,由同一长度的刚性连杆自由连接而成的链的 个非点格算法,其中两个相邻的连杆一道在垂直于点划线轴的平 面内旋转一个在区间〔-Δ卯,△P]中随机选定的角度
事情也可能是这样:不便于(或不可能)在给出的相空间内 对单个自由度均匀抽样。例如,虽然在图24()的左图,从区 间[0,2π]对角度抽样不会遇到什么困难;但在图2.4(e)中,我 们却不能在区问[-∞,∞]内对一变量φ均匀抽样。这个问题出 现于在一个点格上模拟所谓φ模型时: x…+2(2A?+4B)+∑2(,一中 <中;<+0, 其中A,B,是常数,单座点位势v(中)=2A中?+4B中在 A<0,B>0时具有熟悉的双极小的形状。处理这个模型的一种简 单方法,是在超出极小位置(中=±√-A/B)很远的值上 截断中;的容许区间。但是,如果位势很陡峭,这个方法是效率 很低的:大部分时间花在尝试选取试验位形中上,而这些试验 位形却因为跃迁概率太小而被舍弃。如果中,本身是根据一种重 要性抽样方案选取的,就避免了这个问题;这就是说,要建立一 个算法(213,它产生的中的分布正比于概率分布 v(φ p(中)exp(-kaT 最后,图2,4(f)表明,对于随机行走问题,也存在有重要性抽 样方案,以代替迄今讨论的筒单抽样方法。设通过简单抽样建造 出一条长链的一个SAW位形。然后用种种“动力学”算法产生 进一步的位形。例如,在蠕行算法14中,长链一端(选哪 端随机决定)上的终端键被取消,而在另一端的一个随机方向上 加上一个键。当然,这一步试验性变动只有当它不违反SAW约 束时才被执行。 27
这个算法也有一个替代物,如果想要模拟聚合物动力学(例 如,描述一个链在热浴中的弛豫过程的 Rouse模型2.15),这个 替代算法更为现实。这个替代算法允许链上的几组相邻的键局部 重新排列,它们可以随机地跳到点格上的新位置,当然仍然得服 从体斥作用的限制(图2.4(f),中图)(2210。终端键也可以 旋转到新位置上。最后请注意,这些模型都可以有连续的(非点 格的)对应物,例如图2.4(f)右图[217,再者,算法的细节取 决于模拟想达到的目标。例如,如果想研究聚合物熔融的动力 学t217,那么考虑不允许“纠缠”的限制是重要的,即在一条链 的连杆作随机运动时,链必须不与它自身相交,也不与周围别的 链的任何连杆相交。对于使连杆被切掉的运动,这时我们令其跃 迁概率等于零,这样的变动试图永远不会实现。反之,如果我们 主要是对这个模型的静态平衡性质感兴趣,那么最好是这样定义 这场蒙特卡罗模拟的规则,使得能尽快趋于平衡,这意味着应当 不顾不允许发生纠缠的限制,而允许链相交。 显然,要一无遗漏地列举出x→x1这一步的具体内容的全 部各种可能性,以及跃迁概率的详细定义,那是不可能的。蒙特卡 罗方法非常灵活多样,使得它可以应用于多种很不相同的问题。 这显然正是这个方法的主要威力 2.2蒙特卡罗程序的组织和 蒙特卡罗抽样的动力学解释 2.2.1对 Ising模型的模拟的初步讨论 现在,作为一个简单例子,设我们想要实现图2.4(a)的单 个自旋反转 Ising模型的模拟,这怎么做呢? 我们首先必须规定点格的类型和大小,以及必须使用的边界 条件。设我们取一个简单立方点格,大小为LXL×L(即三个 28