由于假定不同的样本i及j是完全独立的，式(3.1-8)中的所有 交叉项消失，可得 =2别-) (3.1-9) =[Kf/u)2)-〈f/w）2] 式(3.1-9)表明I的方差随1/L而变化，不过方差的大小可以由选 择w(x)以使f(x)/e(x)为x的平滑函数而使之大为减小。在理 想情祝下，应使f(x)/w(x)恒定不变，此方差将一同消失。相比之 下，如果w(x)如同在强制性Monte Carlo抽样中一样恒定不变，在 I中的误差则可变得很大。例如，在体积为2的多维构型空间内 抽样，其中只有一小部分∫是可以接受的（例如在前节中，f= 10~260),则在强制性的Monte Carlo抽样中所引起的误差的数量级 为1/(Lf)。由于仅相应于玻尔兹曼常数不为琴的那些构型，式 (3.1-2)中的被积函数才不为签。显然如此抽样才是合理的：进行 非均匀的构型空间抽样，以使得权重函数彻近似地正比于玻尔兹 曼因子。遗憾的是，前面叙述的这种简单重要性抽样方案不能用于 在整个构型空间如式(3：1-2)“对多维积分进行抽样。其理由只是 不知道如何构造由式(3.1-5)到式(3.1-6)的那种变换，使得在 构型空间以正比于玻尔兹曼因子的概率密度来产生各个点。事实 上，对于后一问题的解的一个必要（不过不完全）充分条件是必须 能解析计算狮考察体系的配分函数。当然，如果能对所关心的体系 做到此点，就不会需要计算机模拟了。 3.1.2 Metropolis方法 前一节的最后几行已指明，一般不可能用直接Monte Carlo抽 样计算像ex[-U(r)]dr“这种积分。不过在许多场合，我们 并不对配分函数的构型都分本身，而是对下面的平均值有兴趣： exp[-BU(rN)]A(rN)drN (A) (3.1-10) exp[-BU(rN)]drN 21
因此我们想知道两个积分之比。Metropolis等6]所表明的是有可能 设计一种有效的Monte Carlo方法对这种比值抽样●。为了理解 Metropolis方法，首先让我们更仔细地考察式(3.1-l0)的结构。 此后，将配分函数的构型部分表示为Z Z=jexp[-U(N)]drN (3.1-11) 注意到在式(3.1-10)中比值exP(-βLU)/Z是发现体系处于 附近的概率密度。概率密度可表示为 N(N)=xp[-U(r)】 Z 显然，N(rN)不为零。 假定设法按概率分布N()在构型空间随机产生许多点。这 意味着在点附近单位体积产生的点数为n:,平均来说等于LN (),其中L是所产生的总点数。换句话说 (a)≈是会A(r) (3.1-12) 现在读者大概一定会被式（3.1-12)及3.1,1节中的式(3.1-7) 之间可能产生的差别所图惑。其差别是就式(3.1-7)而言，我们 事先知道在rv附近的假想体积d中抽取一个点的概率，或者是 说exp[-U(rN)]及：2两者都已知。相比之下，在式(3.1 l2)中仅知道exp[-U(N)],即仅知道访问构型空间中不同 点的相对而非绝对概率。这听起来比较抽象，因此让我们惜助于一 简单例于（见图3-1)来阐明此差别。在该图中，我们对比两种测 定尼罗河深度的方法：面积积分法[图3-1(a)]及Metropolis抽 样法，即构造一个重要性权重随机行走[图31(b)]。在常用的 求积方法中，在一组预置点上测定被积函数。由于这些点与被积函 数的值无关，许多点可能位于被积函数值为零的区域。与之不同， 在Metropolis方法中，构造一种随机行走，覆盖被积函数不可忽格 ●Metropolis方法的早期历史的有趣的介绍可在HLAndersen,JStatPhys,43: 731,1986及Wo0d11p.3中查阆。 22
的空间域（即整个尼罗河）。在此随机行走中，凡是带您走出水面 的尝试移动都被拒绝，反之则被接受。经过每一次尝试移动（接受 或拒受)之后，测得水的深度。所有这些测量的（无权重）平均， 获得了尼罗河平均深度的估计值。这就是Metropolis方法的本质。 原则上，常用的求积方法也能给出尼罗河的全部流域的结果。在重 要性抽样方法中却不能得到整个流域的信息，因为此量类似于Z。 尼罗 尼罗河一 (a】 (b) 图3.1测置尼罗河的深度 常提方法(a)和Metropolis方法(b)的比较 下面考虑如何以一个正比于玻尔兹曼因子的相对概率来产生构 型空间的各点。一般的步骤首先在构型空间预备一个体系，称 之为o(老的)，其非零的玻尔兹曼因子为exp[一U(o)]。此构 型可对应于一个无硬核重叠的规则晶体点阵。随之，增加一个小的 随机位移△给o,得到-一个新的尝试构型r,标之以n(新的)。 此新的尝试构型的玻尔兹曼因子为exp[-U(n)]。现必须决定 的是接受或拒受此尝试构型。很多做出此决定的规则都满足这样的 约束，即发现该体系位于构型n的概率正比于N(n)。此处仅讨论 既简单又一般可用的Metropolis方法。 现推导Metropolis方法来确定从构型o至构型n的转移概率 π(o+)。从一个思想实验（实际上是思想模拟）出发较为便利。 我们进行大量的（如M)平行Mionte Carlo模拟，此处M远大于所接 受的构型总数。在任一构型。中的点数示为(o)。我们希望m(o)平 23
均值正比于N(o)。矩阵x(o→)的各元素必须满足一明显的条件：一 旦达到平衡，它们不会破坏平衡。这意味着达到平衡时，导致体系离 开状态o的可接受的尝试移动数目，必须正好等于从所有其他的状态 向状态·的可接受尝试移动数自。也可很便利地附加一更强的条件， 即在平衡时从o至任一其他状态n的可接受移动的平均数目恰好被反 向的移动所抵消。此细致平衡条件意味著着下式 N(o)π(o*n)=N(n)π(n+o) (3.1-13) 转移矩阵π(o→n)的许多可能的形式满足式(3.1-13)。让我们来 看看在实际中π(o→n)是如何构筑的。大家记得一个Monte Carlo 移动由两步组成。首先执行一从状态o至n的尝试移动。用a(o→ n)表示决定执行从i到j的尝试移动的概率转移矩阵，此时a通 常看做是马尔科夫链的基础矩阵[43]。下一步是判定接受或拒受此 尝试移动。今用Pc(o→n)表示接受从o至n的移动的概率， 显然 π(o*n)=a(a*n）XPcc(o+n) (3.1-14) 在最初的Metropolis方法中，a被选为一对称矩阵[Pce(o→n)= Pee(n→o)]。然而在后儿节中，·我们将看到几种a不是对称的例 子。如果a是对称的，可以用Pc(o→n)将式(3.1-13)改写成 N(o)x Pice(o-n)=N(n)x Pacc(n-o) (3.f-15) 由式(3.1-15)得到 Pc（(o→n）=N(nl Pce(n*o） No=exp-U(n以.-U(o)](3.1-16) 同样许多P(on）的不同边择满足此条件[以及明显的条件， 即概率Pe(o→n)不能超过1]。Metropolis等的边择为 Pce(o+n)=N(o)N(n)如果N(n)<N(o) =1 如果N(n)≥N(o)（3.117) 其他的Pac(o→n)的选择也是可行的I9],但Metropolis等的最新 选择较已提出的绝大多数策略更易得到构型空间的有效抽样。 总之，在Metropolis方法中从状态o至状态n的转移概率可 写为 24
x(0+n)=a(o→n) N(n)≥N(o) =a(o→+n)N(n)/N(o)] N(n)<N(o) π(0→0)=1-） ∑π(o→n) (3.1-18) 应予指出，除了矩阵α必须是对称的以外，我们还未指定矩阵a。 这表明在选择各种尝试移动时有相当多的自由度。 至此，仍未解释如何决定一尝试移动是接受或拒受。常用的方 法如下。假定已产生一从状态o至状态n的尝试移动，且U(n)> U(o)。根据式(3.1-16)此尝试移动将以如下概率被接受： Pe(o→n)=exp-β[U(n)-U(o]}<1 为了决定接受或拒受此尝试移动，由区间[0,1]中的均匀分布产 生-随机数，标为R。显然R小于P(o→n)的概率等于Pa (o*n)。如R<Pc(o→n),我们则接受此尝试移动，反之则拒 绝接受。当然，重要的是随机数发生器的确是在区间[0,1]上均 匀地产生随机数。否则Monte Carlo抽样会是偏倚的。对于随机数 的详细讨论可在“数值手册"[30]及Kalos与Whitlock所著的 “Monte Carlo方法"[29]两书中找到。 迄今我们还未提及π(o*n)所必须满足的另一个条件，即各 态历经（即指在有限数目的Monte Carlo步骤中从任一其他点可到 达构型空间的每一可接受点)。显然某些简单的MC方法确保是各 态历经的，但这些方法通常不是最有效的。反之，许多有效的 Monte Carlo方法未被证明是各态历经的，更有甚者，有的被证明 是非各态历经的。解决方案常常是将有效的非各态历经的方法与一 种效率较差但各态历经的方法中的偶然尝试移动相混合。这种方法 总体上则会是各态历经的（至少在原则上）。 3.2基本Monte Carlo算法 很难用抽象的术语来谈论Monte Carlo或分子动力学程序。解 释这类程序如何工作的最好方法是将它们写下来。 绝大多数Monte Carlo或分子动力学程序长度只有几百到几千 25