霍开拓
Monte Carlo 方法在统计物理中的应用 相关书籍和文献 yang元祐
2019-3-13 19:07
阅读:3298

bayes.wustl.edu/Manual/

Equation of State Calculations by Fast Computing Machines

被引数:37469,神作


cftc.cii.fc.ul.pt/ensin

Applications of Monte Carlo methods to statistical physics


assets.cambridge.org/97

A Guide to Monte Carlo Simulations in Statistical Physics, Second Edition


maxwell.ucsc.edu/~peter

Monte Carlo Simulations in Statistical Physics


bookzz.org/book/445395/

A Guide to Monte Carlo Simulations in Statistical Physics


更新


· M. Suzuki, Quantum Monte Carlo Methods in Condensed Matter Physics, World Scientific (Singapore, New Jersey, London, Hong Kong)

· Akira Onuki, Phase Transition Dynamics

Articles:

· Lee, SH., Koo, T.Y. & Jeong, Y.H. Journal of the Korean Physical Society (2016), Ferroelectric domain patterns and domain boundary imprint in the relaxor ferroelectric Pb(Zn1/3Nb2/3)O3-9%PbTiO3 as observed by using piezoresponse force microscopy

· S. C. CHAE, department of physics education, Domain Patterns and Electric Properties at Domain Walls in a Surface Normal to the Direction of Ferroelectric Polarization in h-ErMnO3

· Pereira, A.R. Pires, A.S.T., Gouvêa, M.E. et al. Z. Physik B - Condensed Matter (1992) Dynamical correlation from topological solitons in two-dimensional anisotropic models



蒙特卡洛方法历史

Equation of State Calculations by Fast Computing Machines, N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, Journal of Chemical Physics, 21, 1087-1092 (1953) [PDF]

The Beginning of the Monte Carlo Method, N. Metropolis, Los Alamos Science, 14, 125-130 (1987) [PDF]

Scientific Uses of the MANIAC, H. L. Anderson, Journal of Statistical Physics, 43, 731-748 (1986) [PDF]

蒙特卡洛理论和算法

Progress and Outlook in Monte Carlo Simulations, D. N. Theodorou, Industrial & Engineering Chemistry Research, 49, 3047-3058 (2010) [PDF]

Monte Carlo Simulations, D. J. Earl and M. W. Deem, Methods in Molecular Biology, 443, 25-36 (2008) [PDF]

Monte Carlo Codes, Tools and Algorithms, D. Dubbledam, A. Torres-Knoop and K. S. Walton, Molecular Simulation, 39, 1253-1292 (2013) [PDF]



统计物理学中的蒙特卡罗方法


内容

概述
在统计物理中使用蒙特卡罗方法的一般动机是计算多变量积分。典型的问题始于哈密顿量已知的系统,该系统处于给定温度并遵循玻尔兹曼统计。为了获得某些宏观量的平均值,比如A,一般方法是在所有相空间上计算,为简单起见,使用玻尔兹曼分布计算A的平均值:
{\displaystyle \langle A \rangle = \int _ {PS} A _ {\vec {r}} {\frac {e ^ { - \beta E _ {\vec {r}}}} {Z}} d {\vec {r}}}\\

其中 {\displaystyle E({\vec {r}})= E _ {\vec {r}}} 是由给定状态定义的系统能量。 {\displaystyle {\vec {r}}} 为具有所有自由度的矢量(例如,对于力学系统, {\displaystyle {\vec {r}} = \left({\vec {q}}, {\vec {p}} \right)}  ),   {\displaystyle \beta \equiv 1 / k_ {b} T} 以及 {\displaystyle Z = \int _ {PS} e ^ { - \beta E _ {\vec {r}}} d {\vec {r}}} 是配分函数。

解决这个多变量积分的一种可能方法是精确地列举系统的所有可能构型,并计算平均值。这是在完全可解系统中完成的,并且是在少量粒子构成的简单系统中可模拟的。另一方面,在现实系统中,往往是难以或不可能实现精确枚举的。
对于那些系统,通常采用蒙特卡洛积分。其使用的主要动机是,通过蒙特卡罗积分,误差正比于{\displaystyle 1 / {\sqrt {N}}} ,与积分的维度无关。与蒙特卡洛积分相关的另一个重要概念是重要性采样,它是一种优化计算时间的技术。

重要性抽样


蒙特卡罗积分
{\displaystyle \langle A \rangle = \int _ {PS} A _ {\vec {r}} e ^ { - \beta E _ {\vec {r}}} d {\vec {r}} / Z} 

{\displaystyle \langle A \rangle \simeq {\frac {1} {N}} \sum _ {i = 1} ^ {N} A _ {{\vec {r}} _ {i}} e ^ { - \beta E _ {{\vec {r}} _ {i}}} / Z} 
其中  {\displaystyle {\vec {r}} _ {i}} 从所有相空间(PS)中均匀获取,N是采样点数。
在所有相空间中,一些区域相对于其它区域通常对变量 {\displaystyle A} 的平均值更为重要。特别是那些与剩余的能谱相比足够高的 {\displaystyle e ^ { - \beta E _ {{\vec {r}} _ {i}}}} 值,对积分来说最相关。这样的话,自然要问的是:是否有可能以更密集的方式选择与积分更相关的状态?

那就是使用重要性采样技术
假设  {\displaystyle p({\vec {r}})} 是一种分布,选择与积分更相关的状态.
平均值 {\displaystyle A} 可以改写为
{\displaystyle \langle A \rangle = \int _ {PS} p ^ { - 1}({\vec {r}}){\frac {A _ {\vec {r}}} {p ^ { - 1} ({\vec {r}})}} e ^ { - \beta E _ {\vec {r}}} / Zd {\vec {r}} = \int _ {PS} p ^ { - 1}({ \vec {r}})A _ {\vec {r}} ^ {*} e ^ { - \beta E _ {\vec {r}}} / Zd {\vec {r}}} ,\\ 
其中 {\displaystyle A _ {\vec {r}} ^ {*}} 是考虑重要性概率{\displaystyle p({\vec {r}})}的采样值.

积分可以这样计算
{\displaystyle \langle A \rangle \simeq {\frac {1} {N}} \sum _ {i = 1} ^ {N} p ^ { - 1}({\vec {r}} _ {i} )A _ {{\vec {r}} _ {i}} ^ {*} e ^ { - \beta E _ {{\vec {r}} _ {i}}} / Z}\\  
其中 {\displaystyle {\vec {r}} _ {i}} 使用  {\displaystyle p({\vec {r}})} 分布随机生成。由于大多数时候很难找到一种生成具有给定分布状态的方法,因此必须使用Metropolis算法.


我们知道,x系统的最可能状态是最大化玻尔兹曼分布的状态,一个好的分布  {\displaystyle p({\vec {r}})} ,选择重要性采样是玻尔兹曼分布或经典分布.
{\displaystyle p({\vec {r}})= {\frac {e ^ { - \beta E _ {\vec {r}}}} {Z}}}\\ 
这就是我们要使用的分布。带入之前的求和,
{\displaystyle \langle A \rangle \simeq {\frac {1} {N}} \sum _ {i = 1} ^ {N} A _ {{\vec {r}} _ {i}} ^ {*} }\\ 
要获得给定变量的平均值,使用具有正则分布的 metropolis 算法的过程就是使用metropolis算法来生成由分布 p(r) 给出的状态, 并对 A 实现平均值.

在使用带有正则分布的 metropolis 算法时,必须考虑一个重要问题: 在执行给定的度量(即实现r)时,必须确保该实现与系统之前的状态不相关(否则不会“随机”生成状态). 对于具有能隙的系统,使用正则分布会带来问题,因为系统从先前的状态到脱离关联的新态所需的时间可能趋于无穷大.


Multi-canonical


如前所述,微正则方法有一个大的缺点,这在大多数使用蒙特卡罗积分的系统中都很重要。对于具有“rough energy landscapes”的系统,可以使用Multi-canonical方法。
Multi-canonical方法对重要性采样使用不同的选择:
{\displaystyle p({\vec {r}})= {\frac {1} {\Omega(E _ {\vec {r}})}}}\\ 
其中  {\displaystyle \Omega(E)} 是系统的态密度。这种选择的主要优点是能量直方图是平滑的,即生成的状态在能量上均匀分布. 这意味着,当使用 Metropolis 算法时,模拟不会看到“rough energy landscape”,因为每个能量都被平等对待.
这种选择的主要缺点是,在大多数系统中,  {\displaystyle \Omega(E)} 是未知的. 为了克服这个问题,Wang和Landau算法通常用于在模拟中获取DOS. 请注意,在DOS已知以后,可以为每个温度计算每个变量的平均值,因为状态的生成不依赖于 {\displaystyle \beta} . 

实现
下面,我们看Ising模型。我们考虑一个二维自旋网格,每边都有L个自旋(点阵)。自然有 {\displaystyle N = L ^ {2}} 个自旋,因此,相空间是离散的,其特征是N个自旋.

{\displaystyle {\vec {r}} =(\sigma _ {1},\sigma _ {2},...,\sigma _ {N})} 其中  {\displaystyle \sigma _ {i} \in \{ - 1,1 \}} 是每个格点的自旋. 系统的能量由下式给出 {\displaystyle E({\vec {r}})= \sum _ {i = 1} ^ {N} \sum _ {j \in viz_ {i}}(1-J_ {ij} \sigma _ {i } \sigma _ {j})}\\

其中  {\displaystyle viz_ {i}} 是i的一系列的第一近邻自旋,J是相互作用矩阵(对于铁磁ising模型,J是单位矩阵). 
在这个例子中,目标是获得 {\displaystyle \langle M \rangle}  和 {\displaystyle \langle M ^ {2} \rangle} (例如,为了获得系统的磁化率),因为它可以直接推广到其它可观测量.

根据定义, {\displaystyle M({\vec {r}})= \sum _ {i = 1} ^ {N} \sigma _ {i}}.


首先,必须初始化系统:  {\displaystyle \beta = 1 / k_ {b} T}


有了微观选择,必须采用 metropolis 方法。因为没有正确的方法来选择要选择的那个状态,所以可以尝试在一个时刻翻转一个自旋. 这种选择通常称为单自旋翻转. 执行以下步骤以实现单个测量.


第1步:生成 {\displaystyle p({\vec {r}})}


步骤1.1:执行以下迭代的TT次数:
步骤1.1.1:随机选择一个点阵点(概率为1 / N),称为i,旋转 {\displaystyle \sigma _ {i}}. 
步骤1.1.2:选择一个随机数  {\displaystyle \alpha \in [0,1]}. 
步骤1.1.3:计算试图翻转自旋的能量变化i: {\displaystyle \Delta E = 2 \sigma _ {i} \sum _ {j \in viz_ {i}} \sigma _ {j}}  及其磁化律变化: {\displaystyle \Delta M = -2 \sigma _ {i}} 
步骤1.1.4:如果  {\displaystyle \alpha <\min(1,e ^ { - \beta \Delta E})} ,翻转自旋 ({\displaystyle \sigma _ {i} = - \sigma _ {i}}) ,否则,不翻转。
步骤1.1.5:在自旋翻转的情况下更新宏观量: {\displaystyle E = E + \Delta E}, {\displaystyle M = M + \Delta M}


在T时间之后,系统被认为与其先前状态不相关, 这意味着此时系统处于给定状态的概率遵循玻尔兹曼分布,这是该方法提出的目标.


第2步


步骤2.1:在直方图上保存 M 和 M^2 的值.
最后要注意的是,应该注意到 T 并不容易估计,因为当系统与之前的状态解相关时并不容易。为了克服这一点,人们通常不使用固定的 T,而使用 T 作为隧道时间. 一个隧道时间被定义为步骤1的数量,系统需要从其能量的最小值到其能量的最大值并返回.
在像Ising模型这样的系统中采用自旋翻转选择的这种方法的主要缺点是隧穿时间按照幂律定义为 {\displaystyle N ^ {2 + z}} ,其中z大于0.5,称为临界慢化现象.


适用性


该方法忽略了动力学,这可能是一个主要缺点,或者是一个很大的优点. 实际上,该方法只能应用于静态量,但选择移动的自由度使得该方法非常灵活. 另一个优点是某些系统,例如伊辛模型,缺乏动态描述,仅由能量方法定义; 对于这些,蒙特卡洛方法是唯一可行的方法.


推广


该方法在统计力学中的巨大成功带来了许多的推广, 例如用于优化的模拟退火方法,其中引入了虚拟温度然后逐渐降低.


相关阅读

yang元祐:2维伊辛模型 数值模拟zhuanlan.zhihu.com图标yang元祐:Ising 模型zhuanlan.zhihu.com图标


转载本文请联系原作者获取授权,同时请注明本文来自霍开拓科学网博客。

链接地址:https://wap.sciencenet.cn/blog-566204-1167346.html?mobile=1

收藏

分享到:

当前推荐数:0
推荐到博客首页
网友评论0 条评论
确定删除指定的回复吗?
确定删除本博文吗?