第6章蒙特卡罗方法(3)

1、蒙特卡罗方法（三）6.3 Metropolis 算法算法目的：在积分变量 X 的空间（可能是多维的）内产生一组按概率密度 (X ) 分布的点Metropolis 算法：假想有一个随机行走者在 X 的空间中运动。该随机行走过程相继各步经过的点产生出一个序列：X0, X1, ；随着行走的路程越长，这些点的分布会逼近所要求的分布 (X ) 。注意在这里，行走不是完全无规的。每一步可能达到的位置依赖且仅依赖于上一步的位置。这种随机行走可以用下列跃迁概率描述状态ABCDA1/21/401/4B01/32/30C0100D1/201/20考虑大量行走者从不同的初始点出发在 X 空间独立的随机行走。若 Nn

2、(X ) 是 n 步后这些行走者在 X 点的密度， Nn(Y) 是 n 步后这些行走者在 Y 点的密度，那么在下一步从 X 点走到 Y 点的行走者净数目为X 点上的一个行走者转移到点上的一个行走者转移到Y 的概率的概率当 时，行走者的布局不会有净变化，系统将会达到平衡。这个条件称为细致平衡条件。有多种具体方案给出跃迁概率，其中最重要的是 Metropolis方案细致平衡条件 并不能唯一的定出 跃迁概率为了使行走者的分布达到给定的概率密度分布 (X ) ，需要规定适当的跃迁概率设行走者处于序列中的 Xn 点上，为了产生 Xn+1，行走者迈出试探性的一步到一新点 Xt 。 该新点可以用任何方便的方

3、法选取，如可以在点 Xn 周围的一个边长为 d 的多维立方体中均匀地随机选取。然后按比值来决定是“接受”还是“拒绝”该试探步。Metropolis 的方案Metropolis提出的随机行走的规则这样产生出 X n+1 之后，可以再从X n+1 出发迈出一个试验步，按照同样的过程产生 X n+2 ，不断重复以上步骤，得到整个随机行走序列。如果 r 1，接受该步（即取 X n+1 = X t ）；如果 r 1，则以概率 r 接受这一步：把 r 和一个在0, 1区间上均匀分布的随机数 h 比较，若 h r 就接受这一步（即取 X n+1 = X t ） ，否则就舍弃该步（即取 X n+1 = X n

4、 ）证明 Metropopis方案能导致平衡分布 (X ) 其中 T 是从 X 试探到 Y 的概率，A 是接受这一试探步的概率。根据Metropopis方案，从 X 到 Y 的转移概率为因此，Metropopis随机行走者的平衡分布满足首先试探概率满足若 (X ) (Y )，则 A (YX )=1 ， 细致平衡条件并没有定出跃迁概率的具体形式，因此除了 Metropolis 方案外，还存在其它一些选择，例如 Barker 方法容易证明 Barker 方法确实满足细致平衡条件。其它的选择3. 构成随机行走的点 X0 , X1 , X2 , ，由于产生的方法，彼此不是独立的。1. 随机行走从何出发

5、，即选择何处为 X0 。原则上，任何位置都合适，但实际中，合适的起始点是 值大的地方。Metropolis 方法要注意的地方太小了关联太强，需要很多步才能达到一个统计独立的构型。太大了会使大多数尝试都失败，导致更新缓慢。经验的做法是使得大约的一半试探步入选。2. 步长 d 如何选取？关联函数实际做蒙特卡罗模拟时，需要在生成的点列中进行采样。相邻两个样本点之间应有足够大的的间隔，使得 C(l) 足够小。关联函数定义为自相关函数图例其中 R=(r1, r2, , rN ) 为3N 维坐标矢量，假设粒子间为对势应用实例：单原子气体Metropolis 算法最早被用在经典流体的的模拟中（ Metrop

6、olis, 1953）。这里讨论最简单的经典流体单原子气体。考虑一定体积，温度为T 的单原子气体，这是一个正则系综， 其物理量 A 的期望值为这里没有考虑 U 对速度的依赖，因为速度的分布由麦克斯韦分布给出。 任何依赖于速度的物理量可以直接用该分布解析的计算。下面只考虑依赖于位形的物理量， 如总势能等。 其中Ri 为相空间中依据下面的分布函数抽样得到的样本点用蒙特卡洛方法计算 的方法为其中 hk 为第 k 个方向的步长，、 为 0, 1 区间的随机数很多时候，更新一次构型的只改变一个粒子的坐标效率较高，特别是当系统非常接近平衡态或接近相变时。构型的更新随机的选择第 i 个粒子，其坐标被更新为这

7、步更新被接受的概率为新构型中，只有第 i 个粒子的坐标变化了，因此没有必要为了得到概率 p 而计算整个 (Rn+1)我们将 U 表示为新旧构型的能量差. 构造系统的一个初始态. 随机选择一个粒子 i，产生一个试探步 ri=ri+. 计算这一试探位移引起的能量变化 E. 如果 E 0，生成一个随机数r, 满足0r1. 如果 r exp(- E/KT)，接受这一试探步；否则拒绝. 执行第 (2) 步。 单原子气体的具体算法为了有效的降低计算量，需要对位势引入截断, rij rc 。势函数的截断其截断长度可取为 rc =3，其中 为势函数为 0 时的距离长程相互作用， 例如库仑势，不能被截断，需要额

8、外的方法来处理。例如简单流体中典型的相互作用为 Lennard-Jones 势 步长 h 如果太小： 粒子只能移动很小的距离， 所以需要很多步才能达到一个统计独立的构型。步长 h 如果太大：粒子虽然平均来说可以移动更大的距离，但是绝大多数移动对系统构型的改变过大，以至于能量会显著的增加。 这导致拒绝率过高，从而构型的更新缓慢。步长的选择一个经验法则是接受率平均在 0.4 和 0.6 之间。 对于硬球， 接受率应该低些，约为 0.1。 用来取平均的样本点之间典型的应有 10-15 个间隔样本点。经典的二维 Ising 模型：在一个二维正方格子上，每个格点 i 上有一个自旋，可以取值 +1 或 -

9、1。相邻自旋通过一个交换耦合能 J 相互作用，此外还存在一个外磁场 B。其中 代表对最近邻求和。应用实例Ising 模型系统的哈密顿量为伊辛模型最早被用来研究磁相变，另一个有趣的应用是二元合金。零温下，当没有外磁场时 J 0 ： 如果所有自旋都朝同一方向，系统能量会最低， 对应铁磁态J0：如果相邻两个自旋朝向相反，系统能量会最低。对应反铁磁态交换耦合能 J二维 Ising 模型示意图周期性边条件二维 Ising 模型不同温度下典型的自旋构型二维 Ising 模型的解析结果三维伊辛模型仍然缺乏精确解Lars Onsager 于1944 年得到二维 Ising 模型的精确解，并证明存在相变点，这是统计物理学发展过程中的里程碑。其临界指数为磁化强度磁化率内能比热需要计算的物理量正则系综下磁化强度的平均值被定义为磁化强度可以通过蒙特卡罗方法计算其中 =1，2，M 根据下面的分布函数取样正则系综下物理量的计算其中 为构型 的平均自旋。注意与每一个格点相联系的 可以储存并在模拟中更新。另外，接受概率只有 5 种可能的值，可以事先计算并储存起来，以避免重复进行指数运算。随机的在所有格点生成 1 或 -1， 然后随机的选择一个格点来更新它, 接受率为其中j 对 i 的所有邻居求和。.