3.6在統计力学中的蒙特卡洛方法 假定我们对一个处于熟平衡的恒温(们的体系感兴趣。对该 热力掌问题们敵如下的表述。设有一个包含N个粒子的恒温的 平衡壳系统,我们要计算该系的可观测量A,即该物理量的平 均寶 (4(T)=zJ4G(U( 其中H(x)为系魷的哈密顿量擋述,f()为分布密度函数,Z称为 配分函数,它是归一化常数。 =f((F) 公式中我示在相空间中的庵矢(例如,其坐标为各个数子的空 间位置、动量和自旋普),它给出该状变在相空间中点的坐标。 显蝼上面的公式计算部是涉及很高錐数的积分问题。 在魷计力学的实际问题中,只有京理丸气体、簡谐振子系鴕、 二维 Ising棋型瞢极少数类型的问题可以解析地严格积分求解。 在大多数情况下计算物理量<A>,们只能借助于近似方法求 出。如录用象钟卡洛方法来积分时,只是在把相空间高微化时可 能会引魍误鑾。 如何用蒙卡涪方法來计箕<A>。假宠相应的系是正则系 錄,系统对应于粒子的位形空间參数矢量的哈鲁顿量为 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 模型等极少数类型的问题可以解析地严格积分求解。 在大多数情况下计算物理量< A > ，我们只能借助于近似方法求 出。如果用蒙特卡洛方法来积分时，只是在把相空间离散化时可 能会引起误差。 如何用蒙特卡洛方法来计算< A > 。假定相应的系综是正则系 综,系统对应于粒子的位形空间参数矢量 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大,计算孙到的精度簏而。这个收敛性 是由中心极限定理保证的。由子采用了量耍抽样,我们明显地 提高了数值求解统讣力掌问题<A>的计算效率。 采用 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 越大，计算得到的精度越高。这个收敛性 是由中心极限定理保证的。由于采用了重要抽样法，我们明显地 提高了数值求解统计力学问题< A >的计算效率。 采用 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 常我们对磁化强度的平均值<M(§)>及涨缮<M(S)>-<M(S)>隨系 统的温度和外加叇场的变化感兴慜。它们的讣算公式为 <MS)>z∑MSm0=∑M(sm/∑em, <M(S)>=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 通常我们对磁化强度的平均值< M (S ) > G 及涨落< 2 2 M (S) > − < M (S) > G 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 1 ( ) ( ) ( ) ( ) ( ) ( ) β β β ， ∑ ∑ ∑ − − − − < >= = 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 ( ) ( ) ( ) ( ) ( ) β β β . 按上面公式中的求和来计算的计算量太大，是不可能具体在计算