3.6在统计力学中的蒙特卡洛方法 假定我们对一个处于热平衡的恒温(T)的体系感兴趣。对该 热力学问题我们做如下的表述。设有一个包含N个粒子的恒温的 平衡态系统，我们要计算该系统的可观测量A，即该物理量的平 均值 () ( )( ) xdxHfxAZTA ′′= ∫ ? ′ ? GGG 1 )( . 其中H为系统的哈密顿量描述, ()x′ G ( )xf ′ G 为分布密度函数， Z 称为 配分函数，它是归一化常数。 Z ()()xdxHf ′′= ∫ ? GG . 公式中表示在相空间中的态矢(例如，其坐标为各个粒子的空 间位置、动量和自旋等)，它给出该状态在相空间中点的坐标。 显然上面的公式计算都是涉及很高维数的积分问题。 x′ G 在统计力学的实际问题中，只有象理想气体、简谐振子系统、 二维Ising模型等极少数类型的问题可以解析地严格积分求解。 在大多数情况下计算物理量>< A ，我们只能借助于近似方法求 出。如果用蒙特卡洛方法来积分时，只是在把相空间离散化时可 能会引起误差。 如何用蒙特卡洛方法来计算>< A 。假定相应的系综是正则系 综,系统对应于粒子的位形空间参数矢量x′ G 的哈密顿量为 。如果粒子间的作用力与速度无关，则可以将 中的动能项去掉。这是由于在这时动能项的贡献可以积分积 )(2/)( 1 2 xmpxH N i ii ′Φ+=′ ∑ = GG ()xH ′ G 掉。则在平衡态时其几率分布为波尔兹曼（Boltzmann）分布。 即分布密度函数为 ( ){ } xdZTkxxdxfZxdTxp B GGGGGG 11 )()/(exp))(()(),( ?? Φ?=Φ=， x G 为对动量积分后剩余的相空间坐标，为波尔兹曼常数。上式 中的配分函数为 B k { } ∫ Φ?= )/()(exp TkxxdZ B GG . 从上式中可以看出：所有对应于大能量值的状态x G 对积分贡献都 很小。只有某些状态才贡献很大。因此我们预计在)(x G Φ的平均值 附近分布有很陡峭的峰。采用相空间离散化后的物理量A的系综 平均值表示 () ()() ()() ∑ ∑ = = Φ Φ = n i i n i ii xf xfxA TA 1 1 )( G GG . 我们期望通过随机选择n个状态),...,1( nix i = G ，并对贡献求和的方 法来计算积分。生成的状态越多，物理量A平均值的估计就越精 确。由于相空间是高维的，这就需要产生大量的状态参数，并且 其中大部分的状态对求和的贡献是非常小的。 为了使问题可以有效地进行计算，我们采用重要抽样法的技 术。这种抽样的基本想法是设法产生一个状态的子集合，使其分 布几率为 ( ){ } xdZTkxxdTxp B GGGG 1 )()/(exp),( ? Φ?= . 即取分布概率为系统的热力学平衡态分布。于是系综的物理量A 的平均值就仅仅是对这个状态子集合求平均： ( ∑ = ≈ n i i xA n TA 1 1 )( ) G . n为抽取的状态数。n越大，计算得到的精度越高。这个收敛性 是由中心极限定理保证的。由于采用了重要抽样法，我们明显地 提高了数值求解统计力学问题的计算效率。 >< A 采用Metropolis方法，从任何一个初态出发，进一步生成 一个状态序列。最终生成的状态子集合满足),()( Txpxp GG ≡分布。先 从一个初始状态 0 x G 出发。通过某种抽样方法产生一个状态序列 ???→→→ 321 xxx GGG 。我们规定在单位时间内从系统的一个状态到另一 个状态的过渡几率为 x G x′ G ? ? ? ? ? ? ′ =′ )( )( ,1min),( xp xp xx G w G GG 。抽样方法的选择是至关重 要的。它要能保证抽出的状态子集合满足热力学平衡态分布 )。 ( G xp 从细致平衡条件给出过渡几率只依赖于概率分布的比值这 一事实还可以得到一个重要的结论，即由于状态的分布最终必须 对应于平衡分布( )( )xfZxp GG Φ= ?1 )(，因而比例常数即配分函数Z不会进 入过渡几率。这个结论正反映出这个方法的有用之处。 蒙特卡洛模拟步骤： （a）在相空间中确定一个起始状态 0 x G 。由于马尔科夫链会失去 对初始态的记忆，因而在很大程度上起始状态精确地是什么并不 重要。但是如果初始状态选到与问题无关的那一部分相空间中 时，趋于平衡分布的收敛速度则大大降低。一般选择初始状态处 在分布几率密度最大的区域。 (b)如果已经游走到第n步，现在要游走到第1+n步。产生一个试 探状态或位形 try x G ，使 nn x try x η+=（其中 n η为在间隔 [ ] ?δ δ, 内均匀分 布的随机数）。该状态的选择是：要使δ取得合适。选得太大或 太小，都将很难收敛到平衡分布。选取δ长度的标准是要使1/3 到1/2的试探状态被接受。 (c)计算过渡几率 ( ) tryn xx GG ,w 。 (d)产生一个[0，1]区间的均匀分布随机数 r 。 (e)如果 ( ) tryn xxwr GG ,≤ ，那末接受这一步游走，取 tryn xx GG = +1 。 (f)如果 ( ) tryn xxwr GG ,> ，则把老状态当作新状态，即取 nn xx GG = +1 ，并重新 回到第（b）步。 重复上面的过程，我们就可以完成系统的蒙特卡洛模拟。 这样的模拟过程实际上是对时间的平均。但是, 这是对位形 空间中随时间变化的运动轨道上的物理量所做的平均。统计力学 中的蒙特卡洛模拟方法，按系综不同分为：正则系综蒙特卡洛方 法、微正则系综蒙特卡洛方法、等温等压系综蒙特卡洛方法、巨 正则系综蒙特卡洛方法…。 Ising模型为例来说明正则系综蒙特卡洛方法。 Ising模型是用于解释铁磁性的一个著名的统计格点模型。该 模型的定义如下：令G为一个d 维、共有N个格点的体系； 在每个格子上有一个自旋，可以取朝上或朝下的方向。用自旋 变量来表示 d L= i i s S 如果自旋 . ? ? ? ? = ,1 ,1 i ↓ ↑ 如果自旋 这些自旋之间通过一个交换耦合能J 相互作用。如果还存在一 个外磁场B ，则体系的哈密顿量为 ∑∑∑ == ??= N i i ji j N i i SBSS J 1,1 2 μH . 其中ji,表示只对格点i周围最邻近的格点j求和。μ代表单个自 旋的磁矩。上式中交换耦合能J 是为正时，为铁磁体的模型， 各个自旋倾向于同方向排列；J为负值时，为反铁磁性的模型， 各个自旋倾向于反方向排列。该模型的最大优点就是简单。它忽 略了与格点相关的原子的动能，而仅仅只包括了最相邻原子间的 相互作用能，自旋也仅仅只有两个离散取向。Ising模型尽管简 单，但是利用它仍然可以发现许多有趣的统计性质。 假定相互作用是铁磁性的，即。描述体系性质的配分函 数为 0>J B . ∑ ? = S SH eZ G G )(β 其中)/(1 Tk=β，S }{S= i G 为系统格点上的自旋态位形。任何物理量都 可以由配分函数得到。例如在温度T时的磁化强度为 ∑ ? = ? ? = S SH eSM B Z M G GG )( )( ln1 β β . 其中 M ∑ = = N i i SS 1 )( G 通常我们对磁化强度的平均值< >)(SM G 及涨落< 22 )()( ><?> SMSM GG 随系 统的温度和外加磁场的变化感兴趣。它们的计算公式为 ∑∑∑ ???? =>= S SH S SH S SH eeSMeSMZSM G < G G G G G GGG )()()(1 )()()( βββ ， ∑∑∑ ???? =>=< S SH S SH S SH eeSMeSMZSM G G G G G G GGG )()(2)(212 )()()( βββ . 按上面公式中的求和来计算的计算量太大，是不可能具体在计算 机上计算的。 蒙特卡洛方法则是通过重要抽样，从所有状态的集合中，抽 出一个状态子集合，使得对此子集合中状态的平均与对所有状态 的平均接近，从而算出平均值。然而产生这个状态子集合是通过 一个多次抽样过程来模拟从非平衡态到平衡的弛豫过程来实现 的。 首先随机地给每个格点选取自旋初始值 S ，然后按照顺序, 逐 个地对每个自旋变量通过合适的蒙特卡洛抽样步骤来决定它改 变为另一个状态或者保持不变。 i 对自旋位形抽样的一种基本、常用方法是Metropolis方法。 具体抽样步骤： （1） 选择任意的初始位形S G ={； },...,...,, ssss 21 Ni （2） 按1的等几率，随机抽取一个格点i，将其上的自旋反 向，得到一个新的位形′ N/ S G ={ },...,...,, sss 21 Ni s?； （3） 计算能量差 ( ) ( )SESEE GG ?′=? ，如果0≤?E，则改变有效，取自旋 改变，位形改变 S ii S GG ′ 。这对应于 → )()( SpSp GG >′ 和 W 。 1)( =→ SS ′ GG （4） 如果0>?E，则再产生一个[0，1]区间的随机数r， i （5） 如，则改变仍有效，取自旋改变S E er ?? < β i S GG ′→， （6） 反之（即r），则S E i e ?? ≥ β G 仍保持不变，这对应于 W )}/()(exp{/)}/()(exp{)( TkSHTkSHSS BB GGGG ?′?=′→ . 多次抽样后，就可以逐渐趋于平衡态，得到接近波尔兹曼分布 ( ){ } xdZSHxdSHfZxdTxp G G G G GG 11 )(exp))(()(),( ?? ?== β . 假定我们已经进行了m 次“迭代”，发现系统已经趋于平衡态， 再“迭代”n 次，于是磁化强度的平均值为 ( ∑ + += = nm mi i SM n M 1 1 ) N i ii i . 上面的计算涉及到有2个不同位形的复杂计算，这对于许多大 尺寸的宏观物质的模拟是难于实现的。 为了估计出宏观系统的性质，我们往往给系统强加上周期性 边界条件。即 () ( )LxAxA += GG ， ( ),...,1 di = 其中 L ,为超立方体的线度尺寸。这样就可以 决定下来粒子怎样跨过边界相互作用。在上面的Ising模型中相 互作用仅仅在最临近的格点上的自旋之间，因而这样的周期重复 仅仅是一层格点的复制。周期性边界条件建立起了平移不变性， 这在很大程度上消除了表面效应。 )0,...,0,,0,...,0( L= G ),...,1( diL =