第五章 蒙特卡罗方法(三)ppt课件.ppt
蒙特卡罗方法（三）,5.3 Metropolis 算法,目的：在积分变量 X 的空间（可能是多维的）内产生一组按概率密度 (X ) 分布的点,Metropolis 算法：假想有一个随机行走者在 X 的空间中运动。该随机行走过程相继各步经过的点产生出一个序列：X0, X1, ；随着行走的路程越长，这些点的分布会逼近所要求的分布 (X ) 。,注意在这里，行走不是完全无规的。每一步可能达到的位置依赖且仅依赖于上一步的位置。这种随机行走可以用下列跃迁概率描述,考虑大量行走者从不同的初始点出发在 X 空间独立的随机行走。若 Nn(X ) 是 n 步后这些行走者在 X 点的密度， Nn(Y) 是 n 步后这些行走者在 Y 点的密度，那么在下一步从 X 点走到 Y 点的行走者净数目为,X 点上的一个行走者转移到Y 的概率,当 时，行走者的布局不会有净变化，系统将会达到平衡。这个条件称为细致平衡条件。,有多种具体方案给出跃迁概率，其中最重要的是 Metropolis方案,细致平衡条件 并不能唯一的定出 跃迁概率,为了使行走者的分布达到给定的概率密度分布 (X ) ，需要规定适当的跃迁概率,设行走者处于序列中的 Xn 点上，为了产生 Xn+1，行走者迈出试探性的一步到一新点 Xt 。 该新点可以用任何方便的方法选取，如可以在点 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 ）,证明 Metropopis方案能导致平衡分布 (X ),其中 T 是从 X 试探到 Y 的概率，A 是接受这一试探步的概率。,根据Metropopis方案，从 X 到 Y 的转移概率为,因此，Metropopis随机行走者的平衡分布满足,首先试探概率满足,若 (X ) (Y )， 则 A (XY )=1 ，,在两种情况下，Metropolis 随机行走者的平衡分布都满足,所以行走者最终确实会达到分布 。,若 (X ) (Y )，则 A (YX )=1 ，,细致平衡条件并没有定出跃迁概率的具体形式，因此除了 Metropolis 方案外，还存在其它一些选择，例如 Barker 方法,容易证明 Barker 方法确实满足细致平衡条件。,其它的选择,3. 构成随机行走的点 X0 , X1 , X2 , ，由于产生的方法，彼此不是独立的。,1. 随机行走从何出发，即选择何处为 X0 。原则上，任何位置都合适，但实际中，合适的起始点是 值大的地方。,Metropolis 方法要注意的地方,太小了关联太强，需要很多步才能达到一个统计独立的构型。太大了会使大多数尝试都失败，导致更新缓慢。,经验的做法是使得大约的一半试探步入选。,2. 步长 d 如何选取？,关联函数,实际做蒙特卡罗模拟时，需要在生成的点列中进行采样。相邻两个样本点之间应有足够大的的间隔，使得 C(l) 足够小。,关联函数定义为,自相关函数图例,其中 R=(r1, r2, , rN ) 为3N 维坐标矢量，假设粒子间为对势,应用实例：单原子气体,Metropolis 算法最早被用在经典流体的的模拟中（ Metropolis, 1953）。这里讨论最简单的经典流体单原子气体。,考虑一定体积，温度为T 的单原子气体，这是一个正则系综， 其物理量 A 的期望值为,这里没有考虑 U 对速度的依赖，因为速度的分布由麦克斯韦分布给出。 任何依赖于速度的物理量可以直接用该分布解析的计算。下面只考虑依赖于位形的物理量， 如总势能等。,其中Ri 为相空间中依据下面的分布函数抽样得到的样本点,用蒙特卡洛方法计算 的方法为,其中 hk 为第 k 个方向的步长，、 为 0, 1 区间的随机数,很多时候，更新一次构型的只改变一个粒子的坐标效率较高，特别是当系统非常接近平衡态或接近相变时。,构型的更新,随机的选择第 i 个粒子，其坐标被更新为,这步更新被接受的概率为,新构型中，只有第 i 个粒子的坐标变化了，因此没有必要为了得到概率 p 而计算整个 (Rn+1),我们将 U 表示为新旧构型的能量差,. 构造系统的一个初始态. 随机选择一个粒子 i，产生一个试探步 ri=ri+. 计算这一试探位移引起的能量变化 E. 如果 E 0，生成一个随机数r, 满足0r1. 如果 r exp(- E/KT)，接受这一试探步；否则拒绝. 执行第 (2) 步。,单原子气体的具体算法,为了有效的降低计算量，需要对位势引入截断, rij rc 。,势函数的截断,其截断长度可取为 rc =3，其中 为势函数为 0 时的距离,长程相互作用， 例如库仑势，不能被截断，需要额外的方法来处理。,例如简单流体中典型的相互作用为 Lennard-Jones 势,步长 h 如果太小： 粒子只能移动很小的距离， 所以需要很多步才能达到一个统计独立的构型。步长 h 如果太大：粒子虽然平均来说可以移动更大的距离，但是绝大多数移动对系统构型的改变过大，以至于能量会显著的增加。 这导致拒绝率过高，从而构型的更新缓慢。,步长的选择,一个经验法则是接受率平均在 0.4 和 0.6 之间。 对于硬球， 接受率应该低些，约为 0.1。 用来取平均的样本点之间典型的应有 10-15 个间隔样本点。,经典的二维 Ising 模型：在一个二维正方格子上，每个格点 i 上有一个自旋，可以取值 +1 或 -1。相邻自旋通过一个交换耦合能 J 相互作用，此外还存在一个外磁场 B。,其中 代表对最近邻求和。,应用实例Ising 模型,系统的哈密顿量为,伊辛模型最早被用来研究磁相变，另一个有趣的应用是二元合金。,零温下，当没有外磁场时 J 0 ： 如果所有自旋都朝同一方向，系统能量会最低， 对应铁磁态J0：如果相邻两个自旋朝向相反，系统能量会最低。对应反铁磁态,交换耦合能 J,二维 Ising 模型示意图,周期性边条件,二维 Ising 模型不同温度下典型的自旋构型,二维 Ising 模型的解析结果,三维伊辛模型仍然缺乏精确解,Lars Onsager 于1944 年得到二维 Ising 模型的精确解，并证明存在相变点，这是统计物理学发展过程中的里程碑。其临界指数为,磁化强度,磁化率,内能,比热,需要计算的物理量,正则系综下磁化强度的平均值被定义为,磁化强度可以通过蒙特卡罗方法计算,其中 =1，2，M 根据下面的分布函数取样,正则系综下物理量的计算,其中 为构型 的平均自旋。,注意与每一个格点相联系的 可以储存并在模拟中更新。另外，接受概率只有 5 种可能的值，可以事先计算并储存起来，以避免重复进行指数运算。,随机的在所有格点生成 1 或 -1， 然后随机的选择一个格点来更新它, 接受率为,其中,j 对 i 的所有邻居求和。,. 随机生成一个初始构型. 随机选择一个格点 i. 计算如果将点 i 的自旋翻转引起的能量变化 E. 如果 E 0，生成一个随机数r, 满足0r1. 如果 r exp(- E/KT)，则翻转格点 i 的自旋. 执行第二步。,二维 Ising 模型的具体算法,2020 格子,比热随温度变化曲线,磁化强度随温度变化曲线,
