蒙特卡洛研究伊辛模型小作业例子.pdf
- 文档编号:3431902
- 上传时间:2023-05-05
- 格式:PDF
- 页数:13
- 大小:345.48KB
蒙特卡洛研究伊辛模型小作业例子.pdf
《蒙特卡洛研究伊辛模型小作业例子.pdf》由会员分享,可在线阅读,更多相关《蒙特卡洛研究伊辛模型小作业例子.pdf(13页珍藏版)》请在冰点文库上搜索。
蒙特卡罗方法在二维伊辛模型中的应用姓名:
徐远骥学号:
201328000807008学院:
物理学院指导教师:
周昕2013年11月9日蒙特卡罗方法在二维伊辛模型中的应用一绪论随着计算机技术的不断发展,计算机计算能力的不断提升以及计算机的广泛普及。
计算物理学在科学的探究中发挥着越来越大的作用。
因此我们可以说计算物理学是物理学、数学在过去百余年来取得的巨大成就的基础上,伴随着计算机科学近几十年来突飞猛进的发展而逐渐发展起来的新兴科学分支。
计算物理早已与实验物理和理论物理形成三足鼎立之势,甚至可以说它已成为现代物理大厦的栋梁。
19世纪中叶以前,可以说物理学基本上是一门基于实验的科学。
1862年Maxwell将电磁规律总结为麦克斯韦方程组,进而在理论上语言了电磁波的存在。
这使人们看到了物理理论思维的巨大威力。
从而在理论上开始成为了一门相对独立的物理学分支。
以后到了20世纪初,物理理论经历了两次重大的突破,相继诞生了量子力学和相对论。
理论物理开始成为一门成熟的科学。
物理学研究与计算机和计算机技术紧密结合始于20世纪40年代。
当时正值第二次世界大战时期,美国在研制核武器的工作中,要求准确地计算出与热核爆炸有关的一切数据,迫切需要解决在瞬间内发生的复杂的物理过程的数值计算问题。
然而,采用传统的解析方法求解或手工数值计算时根本办不到的。
这样,计算机在物理学研究中的应用就成为不可避免的事了,计算物理学因此得以产生。
第二次世界大战后,计算机技术的迅速发展又为计算物理学的发展打下了坚实的基础,大大增强了人们从事科学研究的能力,促进了各个学科之间的交叉渗透,使计算物理学得以蓬勃的发展。
计算物理学与理论物理和实验物理有着密切的联系。
计算物理学研究内容涉及物理学的各个领域。
一方面,计算物理学所依据的理论原理和数学方程是由理论物理提供的,其结论还需要实验理论物理来分析检验;另一方面,计算物理学所依赖的数据由实验物理提供的,其结果还要由实验来检验。
对实验物理而言,计算物理学可以帮助解决实验数据的分析、控制仪器设备、自动化数据获取以及模拟实验过程等问题;对理论物理而言,计算物理学可以为理论物理研究提供计算数据,为理论计算提供进行复杂的数值和解析运算的方法和手段。
总之,计算物理学是与理论物理、实验物理互相联系、互相依赖、相辅相成的。
它为理论物理研究开辟了一个新的途径,也对实验物理研究的发展起了巨大的推动作用。
蒙特卡罗方法是一种以概率论和统计理论作为它基础的数值计算方法。
随着计算机技术的发展,这种方法被广泛的应用。
不论是在物理学、数学、计算机科学上都取得了重要的成果。
既然蒙特卡罗是以概率统计理论作为它的理论基础,顾名思义,它所解决的就必然不是像那些具有严格运动形式的时间依赖过程(例如,在经典力学里,如果我们运用牛顿运动方程,就可以精确求解小球在空气中的运动),它所要解决的是那些运动是随机变化的且运动的方式依赖于一系列在模拟中产生的随机数的过程。
对于同样一个初始状态,产生不同序列的随机数我们所得到的结果就不会完全相同,但是这些模拟结果是在一定的统计误差里波动的。
在统计力学中,我们去估计一个模型特定的性质,就需要尝试着在相空间中进行抽样。
很显而易见,由于是蒙特卡罗是一个随机的过程,我们在相空间中移动的路径一定是与模型随时间依赖关系的演化路径是不相同的。
但是在平衡态统计力学中,像具有相互作用的多粒子体系这样的模型,我们的任务是去计算热力学量的平均值。
蒙特卡罗方法就是在考虑合理的统计涨落的情况下可以得到这些模型的相关性质。
所以说,蒙特卡罗的应用范围真的是非常的广泛,很多那些可以用离散化方法近似的模型,蒙特卡罗方法都不失是一种很有效的模拟方法。
二蒙特卡罗方法初步从蒙特卡罗模拟的应用来看,该类型的应用可以分为三种形式:
1.直接蒙特卡罗模拟。
它采用随机数序列来模拟复杂随机过程的效应。
2.蒙特卡罗积分。
这是利用随机数序列计算积分的方法。
积分维数越高,该方法的积分效率就越高。
3.Metropolis蒙特卡罗模拟,这种模拟是以所谓的“马尔科夫”链的形式产生系统的分布序列,该方法可以使我们能够研究经典和量子多体系统的问题。
我们只考虑在经典的情形下相空间C中(C在通常情况下代表很高的维度),在这样高的维度下,蒙特卡罗方法有时候就成为有效求和实际可行的唯一方法。
我们知道配分函数Z是统计问题中的一个重要的函数,在这里我们通常写成在相空间C中具有权重p(x)的积分形式:
=()CZdxpx在经典的系统中,x很显然就表示相空间中的一点。
它的权重常常由玻尔兹曼权重描述为()exp()pxEx,其中()Ex是出于这样位形下的能量。
对于一个物理量A的期望值,我们可以通过遍历相空间中的权重p以及位形中A(x)的平均值来确定,我们写出这样的公式:
1()()pCAdxAxpxZ蒙特卡罗方法的精髓就是将连续问题就行分立化,上式的平均值就可以通过在相空间取一定数量的位形来确定,当我们的位形取得越大的时候,我们的计算结果就越精确,我们就可以得到下面的公式:
11=()MiPMCAAAxM根据中心极限定理的内容,如果我们的位形取得足够的多,那么我们通过蒙卡得到的平均值就会找期望值的附近变化,我们有:
22()()MCPVarAAAAM有时候,我们为了使得蒙特卡罗取样更加合理,就需要选择与()px的分布不同的()x进行计算,这样我们就可以将A写成下面的形式:
()()/()()1()()()/()()(/)/CCCdxAxpxxxAdxAxpxZdxpxxxApp注意,要算上式的平均值,我们就需要对分子、分母分别进行取样计算。
对于上面式子这种最为一般的分布,我们最好是用马科夫(Markov)过程生成位形进行蒙特卡罗的采样。
Markov过程完全由从状态x到状态y这样的转移矩阵xyW来描述。
由几率守恒我们就可以得知1xyyW。
如果满足下面两个条件,我们的Markov过程就将以指数形式从任意一个分部收敛到固定的位形分布()px。
1.各态历经性(Ergodicity):
我们必然可以通过有限次数的Markov过程,从任意一个位形x变换到另一个任意位形y。
也就是说,对于任意的x和y,必然存在着一个整数N,使得nN时,有()0nxyW。
2.细致平衡(detailedbalance)条件:
首先我们要说的是平衡条件。
当系统处于平稳状态时,就要求分布()px满足平衡条件:
()()xyCdxpxWpy,其中()px为转移矩阵xyW的左本征矢。
但是一般情况下我们用的是平衡条件的充分不必要条件进行计算,这就是我们所说的细致平衡条件,形式如下:
()()xyyxWpyWpx蒙特卡罗中的算法中,满足细致平衡条件且最先被提出也是被最广泛应用的算法是Metropolis-Hastings算法。
这种算法的基本思想就是我们认为从位形x到位形y的的更新分别由产生概率propxyW与接受概率accxyW决定。
如果更新被拒绝,那么我们就继续使用原来的位形x。
我们可以写出Metropolis-Hastings算法下的转移矩阵为:
propaccxyxyxyWWW如果我们将接受概率设为:
min1,accxyxyWR其中()()propyxxypropxypyWRpxW这样就可以满足细致平衡条件。
此外,1/yxxyRR。
由此我们就可以得出Metropolis-Hastings算法的一般步骤:
1.首先选取一个试探位置,假定该点位置为trynnxx,其中n为在间隔-,内均匀分布的随机数。
2.计算()()trynfxrfx的数值。
3.如果不等式1r满足,此时()1ntryxx,那就接受这一步游走,并取1ntryxx。
然后返回1,开始对游走到2nx点的试探。
4.如果r1,那么就再另外产生一个0,1区间均匀分布的随机数,如果此时r,那么也还接受这步游走,并取这步游走所到达的点为1ntryxx,然后返回1,开始对游走到2nx点的试探。
如果r,拒绝游走,即仍留在nx位置不变。
5.返回步骤1,重新开始对游走到1nx点的具体位置的又一次试探。
通过上面的步骤,当我们产生了大量的012,.xxx后,才能得到收敛到满足分布f(x)的点集。
三二维伊辛模型的理论介绍在描述磁性系统相变的过程中,有一种最为简单的模型,我们称之为伊辛模型。
例如二维的伊辛模型如图所示,在每个晶格的格点之上都有一个自旋为s的电子,它们自旋的取值即可以向上也可以向下。
对于共有N个格点的晶格系统,若我们只考虑最近邻的相互作用,我们可以写出这个系统的哈密顿量为:
1=JNamijiijiHssHs其中J为耦合系数,s为电子的自旋,可以取1(分别代表自旋向上和自旋向下两个状态),其中代表只对最近邻格点进行求和,为与自旋有关的磁矩,H为磁场强度矢量。
在这里s只代表磁场的大小,而不是量子力学里的算符。
所以伊辛模型在这种情况下实际上是一个准经典的模型。
我们现在只讨论J0的铁磁系统,我们可以写出正则系综的配分函数如下:
12/=.amamNiHkTHkTNssssZee其中Si代表对N个自旋组合的求和。
这样我们就可以求出其它的热力学量,如:
图伊辛模型示意图kTlnZNF2()FETTT=()iiFMsNsH但是由于电子电子之间的耦合作用,要求解这个简单的模型也是十分困难的。
除了一维模型可以严格求解外(1944年昂萨格从平滑的哈密顿量出发对二维伊辛模型进行严格求解)至今人们也不能对三维的情况严格求解。
在这样的情况下,经常我们可以采取某些近似。
我们可以将哈密顿量写成如下的形式:
()amijieffijiJHsHssH其中eff=+hiHH,ijjJhs,其中jjs代表对i的近邻求和。
在这里,我们可以看出,ih与i点附近的电子自旋有关,因此是有涨落效应的。
如果我们做这样的平均场近似,jjss,也就是说将格点i附近的近邻都用一个固定的平均值来代替。
这样,我们就可以写出:
ijjJhs。
我们忽略空间涨落,进一步得到:
jss,=izJhhs其中z代表格点的配位数。
最终,我们就得到了平均场下的哈密顿量的新形式:
1()NamiiHHhs通过平均场近似的方法,我们就将N个自旋耦合的系统近似成了独立的N个电子自旋之和的形式。
由于式中也包含平均自旋,而平均自旋又是我们需要求的,这样就构成了一个自洽求解的方程。
我们参阅相关的统计物理学的书籍,可以知道确定s的自洽方程为:
tanh()HzJsskTkT四蒙特卡罗方法模拟二维伊辛模型在这里我们采用蒙特卡罗方法对二维伊辛模型的顺磁铁磁转变就行模拟。
我们取约化能量单位为0/1BJkT。
这样我们可以将计算公式简写为如下的形式。
总能量为:
()TUTE,总磁矩为:
()=iTiMT,总能量的方差为:
22var()TTEEE,总磁矩的方差为:
22var()TTMMM。
我们取了100*100的二维平方格点,将温度T从1每隔0.1一直取到3.5,观察总能量,总磁矩与比热随温度T的变化关系。
(注意:
在这里T并不是真的温度,而是约化单位中T0的倍数)。
首先我们要讨论,什么时候我们的取样才不依赖于初始温度,即在第一步的热平衡中,我们的预热是否充分,在我的模拟中,我们热平衡的步数取为一百万步(即程序中的thermal=1000000),我们观察从初始到这一百万个蒙特卡罗步中,单个格点上的能量与磁矩随蒙特尔卡罗步数的变化趋势,就可以大体知道我们的取样是否充分。
如上图所示,从单个格点上的能量来看,在20万步左右的情况下,总能量就在一个很小的范围内波动。
在从下面的磁矩上看,我们可以看到,开始温度较低的时候如T=1时,在20万步的情况下,系统就进入了平衡状态,单个格点上的磁矩在1附近几乎不发生变动。
在T升高的过程中,单个格点上的磁矩越来越小,对应的总磁矩为趋近与零的情形。
通过上的分析来看,从总体来说,我们取的100万步的蒙特卡罗步是一个比较合理的值。
下面我们分析平均接收概率随温度的变化情况,从下图我们可以看出,在温度很低的时候,接受率非常的低,如T=1时,接受率仅为0.0007243。
对于这样的模拟,要达到我们设定的足够多的样本数,就需要花费非常多的时间。
在本次模拟中,我们采取这样的策略,我们固定要抽样的样本数,而不管接受率的大小,通过延长模拟时间来得到可靠的样本。
在我们的模拟中,在接受率较大的时候如T=3时,我们进行10万步的抽样,所花时间也就在10分钟左右,然而在T=1时,我们要花费数小时的时间。
下面,我们将所得到的的不同T值时的数据进行分析,我们发现在温度较低时,系统处于铁磁状态,随着温度的升高,系统的磁化率越来越低,最后处于顺磁状态。
这T=2.25附近很明显存在一个相变的过程,我们称这样的相变为铁磁顺磁相变。
这种低温下的不需要外接作用的磁化就称之为自发磁化。
顺磁-铁磁转变的温度就为临界温度。
这这样的磁化下,其实总自旋是有正负之分的,分别为自旋向上和自旋向下,在本程序中统计时我们取的是绝对值,强制认为向上和向下是空间对称的。
其实上,这种向上向下的区分产生真是一种自发的对称性破缺。
根据朗道二级相变的平均场理论,二级相变正是存在这样的对称性破缺。
我们也知道,对于二级相变,在相变点,两相的化学势和化学势的一级偏微商全部相等,但化学势的二级偏微商不相等,这样我们的定压热容量就会在相变点存在一个发散的行为,如图所示:
同时昂萨格1944年通过理论求出了二维伊辛模型的严格解,他指出,对二维正方形点阵,转变温度为:
2.269CJTk。
我们模拟中取J=2,k=1。
可以看出我们的转变温度与理论值是非常一致的。
通过多次模拟我们发现,在温度较高的情况下,我们的结果通常是可信的,但是温度较低时,我们有时很难达到一个平衡状态,100万步是不够的,因为在低温T很小时,接收概率小,虽然我们也取了同样多的样本,但由于位形变化的幅度很小,通过多次变化得出的位形其实没有太大的变化,因此位形之间的关联是很大的,并不能得到足够多独立的样本。
在这样的情况下,我们就需要增加两次抽样中所需要间隔的位形数目,如下图所示,我们就可以清楚的看到,在低温下,如果抽取的样本被trap在几个固定的位形时,这个抽样样本之间很强的关联就会导致所得到的结果非常糟糕。
如T=1.1与T=2附近的那几个值都是不可用的值。
那么如何证明我们的说法是否正确呢,下面我们将热平衡时的自旋平均值与正常情况下的自旋平均值做一个比较。
如左下图所示,正常情况下,确实我们可以看出,由于样本间的关联很小,在适当的步数之后,我们可以达到一个真正的热平衡状态。
这就是图中红色线所表示的平均单个格点上的自旋磁矩与蒙特卡罗步的关系。
在40万步时我们已经进入了平衡状态。
而在样本关联很大的情况下,我们的位形就被束缚在了一个位形的附近,我们通过分析抽样的样本发现,在过了很长的蒙特卡罗步之后,我们的样本才开始进入一个热平衡的状态,这就是图中绿色线条所表示的情况。
通过这样的错误样本我们当然会得到很差的结果。
五总结在本次模拟中,我们用蒙特卡罗方法对二维平方格点的伊辛模型进行模拟。
我们将温度T从1每隔0.1取到3.5。
我们分析了所有的样本,看到了二维伊辛模型的顺磁铁磁转变,并从理论上分析了二维伊辛模型的一些性质。
最终我们发现我们的模拟结果符合理论上的预期。
当然我们也对模拟中产生的问题进行了讨论。
如低温下接受概率很低的问题,如何判断系统进入平衡态的问题,以及对于关联很强的样本导致结果不准确的问题。
通过这次模拟,对二维伊辛模型理论,我有了更加深刻的理解,也充分认识到了计算机模拟在现代物理问题解决中的重要性。
更加重要的是我们要从大量的数据中分析问题,找到有用的结果,并且能对我们的结果进行一个合理的估计。
六参考文献1马文淦.计算物理学.北京:
科学出版社M.2012:
53-54.2林宗涵.热力学与统计物理学.北京:
北京大学出版社M.2007:
172-178.3谭浩强.C+程序设计.北京:
清华大学出版社M.2010.4DaanFrenkel&BerendSmit.MolecularSimulation:
fromalgorithmtoapplications.SanDiegoSanFranciscoNewYork:
AcademicPressM.2002.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 蒙特卡洛 研究 模型 作业 例子
![提示](https://static.bingdoc.com/images/bang_tan.gif)