3.6在統计力学中的蒙特卡洛方法 假定我们对一个处于熟平衡的恒温(们的体系感兴趣。对该 热力掌问题们敵如下的表述。设有一个包含N个粒子的恒温的 平衡壳系统,我们要计算该系的可观测量A,即该物理量的平 均寶 (4(T)=zJ4G(U( 其中H(x)为系魷的哈密顿量擋述,f()为分布密度函数,Z称为 配分函数,它是归一化常数。 =f((F) 公式中我示在相空间中的庵矢(例如,其坐标为各个数子的空 间位置、动量和自旋普),它给出该状变在相空间中点的坐标。 显蝼上面的公式计算部是涉及很高錐数的积分问题。 在魷计力学的实际问题中,只有京理丸气体、簡谐振子系鴕、 二维 Ising棋型瞢极少数类型的问题可以解析地严格积分求解。 在大多数情况下计算物理量,们只能借助于近似方法求 出。如录用象钟卡洛方法来积分时,只是在把相空间高微化时可 能会引魍误鑾。 如何用蒙卡涪方法來计箕。假宠相应的系是正则系 錄,系统对应于粒子的位形空间參数矢量的哈鲁顿量为 H(x)=∑n212m+0(x)。如录粒子间的作用力与度元关,则可以将 H(x)中的动能项去掉。这是由于在这时动能项的贡就可以积分积
3.6 在统计力学中的蒙特卡洛方法 假定我们对一个处于热平衡的恒温(T)的体系感兴趣。对该 热力学问题我们做如下的表述。设有一个包含 N 个粒子的恒温的 平衡态系统，我们要计算该系统的可观测量 A，即该物理量的平 均值 A T = Z A( ) x′ f (H(x′))dx ∫ − ′ Ω 1 G G G ( ) . 其中H(x′)为系统的哈密顿量描述, G f (x′) G 为分布密度函数， Z 称为 配分函数，它是归一化常数。 Z = f (H( ) x′ )dx′ ∫ Ω G G . 公式中 表示在相空间中的态矢(例如，其坐标为各个粒子的空 间位置、动量和自旋等)，它给出该状态在相空间中点的坐标。 显然上面的公式计算都是涉及很高维数的积分问题。 x′ G 在统计力学的实际问题中，只有象理想气体、简谐振子系统、 二维 Ising 模型等极少数类型的问题可以解析地严格积分求解。 在大多数情况下计算物理量 ，我们只能借助于近似方法求 出。如果用蒙特卡洛方法来积分时，只是在把相空间离散化时可 能会引起误差。 如何用蒙特卡洛方法来计算 。假定相应的系综是正则系 综,系统对应于粒子的位形空间参数矢量 x′ G 的哈密顿量为 。如果粒子间的作用力与速度无关，则可以将 中的动能项去掉。这是由于在这时动能项的贡献可以积分积 ( ) / 2 ( ) 1 2 H x p m x N i i i ′ = ∑ + Φ ′ = G G H(x′) G
掉。则在平衡恋时其几亭分布为波尔兹曼( Boltzmann)分布。 即吩布寧度函数为 p(,T)d=(Z)f((x)d=exp()k)Z)4, 为对动量积分后剩汆的相空间坐标,k为波尔兹受常数。上式 中的配分函教为 exp(Φ(x)(kBT)}. 从上式中可以看出:所有对疝于大能量值的状对积分贡觥部 很小。只有某些状亮攻很大。因此我们预计在Φ(x)的平均值 附近分布有很陡峭的峰。采用相空间高化后的物理量A的系综 平均值表示 ∑4()(() AT ∑f(G) 我们期望通过随机选择n个状x(i=1…,n),外对贡就求和的方 法來计算积分。生成的状虍多,物理量A平均值的佑计就 确。由于相空间是商维的,这就鼎产生大量的状參数,并且 其中大部分的状虍对求和的族是非常小的。 为了良问题可以有效地进行计算,我们采用置要抽样法的技 米。这种抽祥的基本法是设法产生一个状充的子合,使其分 布几率为 p(x,T)d=exp(ap()/(kgT)XZ)dx 即取分布概率为系统的熟力学平衡庵分布。于是系缭的物理量A 的平均寶就仅仅是对这个状子合求平均:
掉。则在平衡态时其几率分布为波尔兹曼（Boltzmann）分布。 即分布密度函数为 p x T dx Z f x dx { (x) k T } Z dx B G G 1 G G G 1 G ( , ) ( ) ( ( )) exp /( ) ( ) − − = Φ = − Φ ， x G 为对动量积分后剩余的相空间坐标， 为波尔兹曼常数。上式 中的配分函数为 B k { } Z = ∫ dx exp − Φ(x) /(k T ) B G G . 从上式中可以看出：所有对应于大能量值的状态x G 对积分贡献都 很小。只有某些状态才贡献很大。因此我们预计在 (x) G Φ 的平均值 附近分布有很陡峭的峰。采用相空间离散化后的物理量 A 的系综 平均值表示 ( ) ( ) ( ) ∑ ( ) ( ) ∑ = = Φ Φ = n i i n i i i f x A x f x A T 1 1 ( ) G G G . 我们期望通过随机选择 n 个状态 x (i 1,..., n) i = G ，并对贡献求和的方 法来计算积分。生成的状态越多，物理量 A 平均值的估计就越精 确。由于相空间是高维的，这就需要产生大量的状态参数，并且 其中大部分的状态对求和的贡献是非常小的。 为了使问题可以有效地进行计算，我们采用重要抽样法的技 术。这种抽样的基本想法是设法产生一个状态的子集合，使其分 布几率为 p x T dx { (x) k T } Z dx B G G G 1 G ( , ) exp /( ) ( ) − = − Φ . 即取分布概率为系统的热力学平衡态分布。于是系综的物理量 A 的平均值就仅仅是对这个状态子集合求平均：
4)=∑() n为抽取的状数。n大,计算孙到的精度簏而。这个收敛性 是由中心极限定理保证的。由子采用了量耍抽样,我们明显地 提高了数值求解统讣力掌问题的计算效率。 采用 Metropolis方法,从任何一个初庵出墩,进一步生成 个状序列。最终生成的状壳子集合演足p(x)=p(,7)分布。先 从一个初始状忞出敬。過过神抽祥方法产生一个状序列 →2→1→…。我们规定在单位时间内从系統的一个状走到另 个状庵氵的过渡几率为。,)2=mL3l。抽样方法的选择是至关量 要的。它要能保证抽出的状态子集合演足熟力学平衡走分布p(x)o 从细政予衡条件給出过渡几率只依赖于概率分布的比值这 一喜还可以得到一个量耍的结论,即由于状壳的分布最终必须 对应于平衡分布p()=zf((),因而比例常数即配分函数z不拿进 尅渡几率。这个结论正凤映出逭个方法的有用之处。 特卡洛拟步骤 (a)在相空间中确定一个始状走。由于马尔科失链会失去 对物始的记忆,因而在很大程度上兔始状变崭确地是什么外不 县。但是如果初始状选到与问逦无关的那一部分相空间中 时,于平衡分布的收敛遠度则大大降氐。一般选择初始状充处 在分布几率譽度最大的区域。 (b)如果已经游走到第n步,现在要游先到第n+1步。产生一个试 探状变成位形,良x=xn+n(其中n为在间隔[-。内均訇分
∑ ( = ≈ n i i A x n A T 1 1 ( ) ) G . n 为抽取的状态数。n 越大，计算得到的精度越高。这个收敛性 是由中心极限定理保证的。由于采用了重要抽样法，我们明显地 提高了数值求解统计力学问题的计算效率。 采用 Metropolis 方法，从任何一个初态出发，进一步生成 一个状态序列。最终生成的状态子集合满足 p(x) p(x,T) G G ≡ 分布。先 从一个初始状态 0 x G 出发。通过某种抽样方法产生一个状态序列 → → → ⋅⋅⋅ 1 2 3 x x x G G G 。我们规定在单位时间内从系统的一个状态 到另一 个状态 的过渡几率为 x G x′ G ′ ′ = ( ) ( ) ( , ) min 1, p x p x w x x G G G G 。抽样方法的选择是至关重 要的。它要能保证抽出的状态子集合满足热力学平衡态分布 ( )。 G p x 从细致平衡条件给出过渡几率只依赖于概率分布的比值这 一事实还可以得到一个重要的结论，即由于状态的分布最终必须 对应于平衡分布 p x Z f ( (x)) G G = Φ −1 ( ) ，因而比例常数即配分函数Z 不会进 入过渡几率。这个结论正反映出这个方法的有用之处。 蒙特卡洛模拟步骤： （a）在相空间中确定一个起始状态 0 x G 。由于马尔科夫链会失去 对初始态的记忆，因而在很大程度上起始状态精确地是什么并不 重要。但是如果初始状态选到与问题无关的那一部分相空间中 时，趋于平衡分布的收敛速度则大大降低。一般选择初始状态处 在分布几率密度最大的区域。 (b)如果已经游走到第n步，现在要游走到第n +1步。产生一个试 探状态或位形 try x G ，使 n n x try x = +η （其中ηn为在间隔[−δ ,δ ]内均匀分
布的隨机数)。该状恋的选择是:要俍δ取得合遁。选憬太大或 太小,将很难收敛到平衡分布。选取δ长度的橛准是要1/3 到1/2的谜捍状态被接曼。 (c)计算过渡几率vnx (d)产生一个[0,1]区间的垍訇分布隨机数r。 (e)如果,sm,元),那來换受这一步游走,取x=元。 (f)如果,>M(,x),则把老状庵当作新状走,即取=,外堂新 回到第(b)乡。 复上面的过程,就们就可以兜成系统的毀特卡洛樓拟。 这禅的模拟过程实际上是对时间的平均。卫是,这是对位形 空间中时间变化的返动軌熴上的物理量所僦的平均。筑讣力 中的特卡涪模拟方法,换系錄不同分为:正则系錄毀特卡洛方 法、微正则系缭靈特卡洛方法、普温普压系綠象兮卡洛方法、亘 正则系缭象特卡洛方法 Ising棋型为例来说明正则系鲦蒙特卡洛方法。 Ising模型是用于解尋铁礅性的一个着名的鲩计格点棋型。该 模型的定义如下:令G=为一个d维、共有N个格点的体系; 在每个格子i上有一个自旋,可以取朝上成眀下的方向。用自旎 变量s来表示 S 1,如果自旋↑ 1,如果自旋↓ 这些自旎之间过一个交换耦合能丁相互作用。如最还存在
布的随机数）。该状态的选择是：要使δ 取得合适。选得太大或 太小，都将很难收敛到平衡分布。选取δ 长度的标准是要使 1/3 到 1/2 的试探状态被接受。 (c)计算过渡几率 ( ) n try x x G G w , 。 (d)产生一个[0，1]区间的均匀分布随机数 r 。 (e)如果 ( ) n try r w x x G G ≤ , ，那末接受这一步游走，取 n try x x G G +1 = 。 (f)如果 ( ) n try r w x x G G > , ，则把老状态当作新状态，即取 n n x x G G +1 = ，并重新 回到第（b）步。 重复上面的过程，我们就可以完成系统的蒙特卡洛模拟。 这样的模拟过程实际上是对时间的平均。但是, 这是对位形 空间中随时间变化的运动轨道上的物理量所做的平均。统计力学 中的蒙特卡洛模拟方法，按系综不同分为：正则系综蒙特卡洛方 法、微正则系综蒙特卡洛方法、等温等压系综蒙特卡洛方法、巨 正则系综蒙特卡洛方法…。 Ising 模型为例来说明正则系综蒙特卡洛方法。 Ising 模型是用于解释铁磁性的一个著名的统计格点模型。该 模型的定义如下：令G 为一个 d 维、共有 N 个格点的体系； 在每个格子 上有一个自旋，可以取朝上或朝下的方向。用自旋 变量 来表示 d = L i i s S 如果自旋 . − = 1, 1, i ↓ ↑ 如果自旋 这些自旋之间通过一个交换耦合能 J 相互作用。如果还存在一
个外礅场B,则体系的哈密顿量为 H=-7∑S∑S-B∑S 其中示只对格点周图最邻近的格点j求和。代表单个自 旋的磁矩。上式中交换耦合能J是为正时,为铁礅体的模型, 各个自旋倾向于同方向排列;丁为负值时,为反铁礅性的模型, 各个自旋傾向于反方向排列。该棋的最大优痕就是简单。它忽 略了与格点相头的原子的动能,而仅仅只包括了量相邻原子间的 相互作用能,自旋也仅仅只有两个离散取向。 Ising模型恳蒿 ,但是利用它仍然可以墩现许多有越的統计性质。 假定相工作用是铁礅性的,即J>0。描述体系性孜的配分函 教为 z=e-ph( 其中B=1k2D),§={S}为系格点上的自旋壳位形。任何物理量都 可以由配分函教得到。例如在温度T时的叇化强度为 M=1an2=∑M(sm 其中 (S)=∑S 常我们对磁化强度的平均值及涨缮-隨系 统的温度和外加叇场的变化感兴慜。它们的讣算公式为 z∑MSm0=∑M(sm/∑em, =z-EM(S)e H(S)=2M"(S)e-pH(5)/2e"(5) 驶上窗公式中的求和来计犷的计量太大。是不可能具体在计
个外磁场 B ，则体系的哈密顿量为 ∑ ∑ ∑ = = = − − N i i i j j N i i S S B S J 2 1 , 1 H µ . 其中 i, j 表示只对格点i周围最邻近的格点 j 求和。µ 代表单个自 旋的磁矩。上式中交换耦合能 J 是为正时，为铁磁体的模型， 各个自旋倾向于同方向排列；J 为负值时，为反铁磁性的模型， 各个自旋倾向于反方向排列。该模型的最大优点就是简单。它忽 略了与格点相关的原子的动能，而仅仅只包括了最相邻原子间的 相互作用能，自旋也仅仅只有两个离散取向。Ising 模型尽管简 单，但是利用它仍然可以发现许多有趣的统计性质。 假定相互作用是铁磁性的，即 。描述体系性质的配分函 数为 J > 0 B ∑ . − = S H S Z e G G β ( ) 其中β = 1/(k T)，S ={S }i G 为系统格点上的自旋态位形。任何物理量都 可以由配分函数得到。例如在温度 T 时的磁化强度为 ∑ − = ∂ ∂ = S H S M S e B Z M G G G ( ) ( ) 1 ln β β . 其中 M ∑= = N i i S S 1 ( ) G 通常我们对磁化强度的平均值 G 及涨落 − G G 随系 统的温度和外加磁场的变化感兴趣。它们的计算公式为 ∑ ∑ ∑ − − − − >= = S H S S H S S H S M S Z M S e M S e e G = = S H S S H S S H S M S Z M S e M S e e G G G G G G G G G 2 1 2 ( ) 2 ( ) ( ) ( ) ( ) ( ) β β β . 按上面公式中的求和来计算的计算量太大，是不可能具体在计算
机上计算的。 象特卡洛方法则是逦过量要抽禅。从所有状庵的集合中。抽 出一个状态子合,得对此子集合中状亮的平均与对所有状亮 的平均换近,从而算出平均。粪而产生这个状子集合是遭过 一个多次抽样程来棋拟从非平衠恋到平衡的弛豫过程来实觊 的 首先随机地給每个格京选取自捩初始s,然后换照顺序,遝 个地对每个自旋变量过合适的象特卡洛抽禅步骤来决定它欧 变为另一个状成者保不变。 对自旋位形抽禅的一种基本、常用方法是 Metropolis方法。 具体抽祥步骤 (1)选择任宽的初始位形={s,s2,s,5}; (2)换1N的等几亭,机抽取一个格点i,将其上的自旋反 向,得到一个新的位形={s1,52-s灬s} (3)计算能量△E=s)-E(,如果AE≤,则欧变有效,取自旎 改变,位形改变→。这对应于p(S)>p(5)和ws→s (4)如果△E>0,则再产生一个[0,1区间的随机数r, (5)如nS)=exp(H(S)/(kgT))/exp-H(S)/(kBT)) 多次抽样后,就可以還漸趋于平衡尧,得到接近波尔兹曼分布 p(,)=(2)f(H(5)=pmB)
机上计算的。 蒙特卡洛方法则是通过重要抽样，从所有状态的集合中，抽 出一个状态子集合，使得对此子集合中状态的平均与对所有状态 的平均接近，从而算出平均值。然而产生这个状态子集合是通过 一个多次抽样过程来模拟从非平衡态到平衡的弛豫过程来实现 的。 首先随机地给每个格点选取自旋初始值S ，然后按照顺序, 逐 个地对每个自旋变量通过合适的蒙特卡洛抽样步骤来决定它改 变为另一个状态或者保持不变。 i 对自旋位形抽样的一种基本、常用方法是 Metropolis 方法。 具体抽样步骤： （1） 选择任意的初始位形S G ={s1 ,s2 ,...,si ,...sN }； （2） 按1 的等几率，随机抽取一个格点i，将其上的自旋反 向，得到一个新的位形 ′ / N S G ={s ,s ,..., s ,... } 1 2 i N − s ； （3） 计算能量差 E E(S ) E(S) G G ∆ = ′ − ，如果∆E ≤ 0，则改变有效，取自旋 改变，位形改变Si i S G G → ′。这对应于 p(S ) p(S ) G G ′ > 和W (S → S ′) = 1。 G G （4） 如果∆E > 0，则再产生一个[0，1]区间的随机数ri， （5） 如 ，则改变仍有效，取自旋改变S E r e − ∆ < β i S G G → ′， （6） 反之（即r ），则S E i e − ∆ ≥ β G 仍保持不变，这对应于 W (S S ) exp{ H (S ) /(k T )} / exp{ H (S ) /(k T )} B B G G G G → ′ = − ′ − . 多次抽样后，就可以逐渐趋于平衡态，得到接近波尔兹曼分布 p x T dx Z f H S dx { H(S )} Z dx G G G G G G 1 1 ( , ) ( ) ( ( )) exp ( ) − − = = − β
假定我们已经进行了Ⅲ次“选代”,发现系統已经于平衡, 再“迭代》n次。于是磁化强度的平均值为 上面的计算涉及到有2个不同位形的复杂计算,这对于许多大 尺寸的宏观物孜的棋拟是难于实现的。 为了佑计出宏观系統的性放,我仉往往给系統强加上周期性 逸界条件。即 AG)=A(x+L),(i=1,,d) 其中L=(0,、0L10.0),L(=1,d为繩立方体的线度尺寸。这禅就可以 决定下来粒于夜禅跨过边界相互作用。在上面的 Ising桃型中相 互作用仅仅在量临近的格底上的自旋之间,因而这样的周期量复 仅仅是一层格京的复侧。周期性边界永件~文兔丁平移不变位, 这在很大程度上消除了表面效应
假定我们已经进行了 m 次“迭代”，发现系统已经趋于平衡态， 再“迭代”n 次，于是磁化强度的平均值为 ∑ ( + = + = m n i m M Si n M 1 1 ) N i i i i . 上面的计算涉及到有2 个不同位形的复杂计算，这对于许多大 尺寸的宏观物质的模拟是难于实现的。 为了估计出宏观系统的性质，我们往往给系统强加上周期性 边界条件。即 A( ) x = A(x + L ) G G ， (i = 1,..., d ) 其中L , 为超立方体的线度尺寸。这样就可以 决定下来粒子怎样跨过边界相互作用。在上面的 Ising 模型中相 互作用仅仅在最临近的格点上的自旋之间，因而这样的周期重复 仅仅是一层格点的复制。周期性边界条件建立起了平移不变性， 这在很大程度上消除了表面效应。 = (0,...,0, L ,0,...,0) G L (i = 1,..., d)