第五章-蒙特卡罗方法(三)

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