蒙特卡诺模拟设计内附详细程序设计例子.docx

1、蒙特卡诺模拟设计内附详细程序设计例子蒙特卡诺（Monte-Carlo）方法（一）Monte-Carlo模拟简介 Monte-Carlo方法亦称为随机模拟（Random Simulation）方法，有时也称作随机抽样（Random Sampling）技术或统计试验（Statistical Testing）方法。它的基本思想是，为了求解数学、物理、工程技术以及生产管理等方面的问题，首先建立一个概率模型或过程的观察或抽样试验来计算所求参数的统计特征，最后给出所求解的近似值。是由Metropolis在二次世界大战期间提出的Manhattan计划，用于研究与原子弹有关的中子输运过程中提出的。应用Mont

2、e-Carlo方法解决实际问题：1、对求解的问题建立简单而又便于实现的概率统计模型，使所求的解恰好是所建立模型的概率分布或数学期望。2、根据概率统计模型的特点和计算实践的需要，尽量改进模型，以便减小方差和降低费用，提高计算效率。3、建立对随机变量的抽样方法，其中包括建立产生伪随机数的方法和建立对所遇到的分布产生随机变量的随机抽样方法。4、给出获得所求解的统计估计值及其方差或标准误差的方法。Monte-Carlo模拟在物理研究中的作用：统计物理学与具有多自由度的系统打交道。统计物理学的一个任务是从模型的哈密顿量计算出所要的各种平均性质，例如每个自由度的平均能量E或平均磁化强度M， EtN, M=

3、tN (2.7.1)其中t代表热平均。A(x)表示任何可观察量，例如A=H,等等，x是相空间中的矢量，代表所考虑的自由度的一组变量的集合，在我们的研究范畴内，x=(),表示系统的自旋组态。对下式式x=()的热平均在正则系综中由下式定义： T =, (2.7.2) Z=. 归一化的Boltzmann因子 P(x)= (2.7.3)起着一个概率密度的作用，它描写位形x出现在热平衡中的权重。MonteCarlo模拟实际上是一个重要性抽样，重要性抽样的想法是从对一个量（例如磁化强度）有最大贡献的那些场合进行抽样。现在我们考虑一个抽样过程，在这个过程中，我们根据某一个概率P(xl)来选取相空间点xl。然

4、后取这一组xl来求热平均， (2.7.5)式中P(xl)的一种简单而且最自然的取法是 ；这时Boltzmann因子完全相消，于是(1)式便化为简单算术平均： (2.7.6)为了找到一种抽样方法，来实际实现这个所谓重要性抽样，Metropolis等人提出了下述想法：可以建造一个Markov过程，而不要彼此独立无关地选取相继诸状态xl，过程中通过一个适当的跃迁概率W(xlxl+1)，每一个状态xl+1由前一个状态xl得到。同时，他们还指出，有可能这样选取跃迁概率W，使得在M的极限下，Markov过程产生的状态分布函数P(xl)趋于所要的平衡分布 (2.7.7)而达到这样一个结果的充分条件是加上细致

5、平衡条件 (2.7.8)上式意味着，变迁的跃迁概率和反方向变迁的跃迁概率之比，只会依赖于能量的变化 (2.7.9)由此可见，原来平衡分布时的概率密度函数的分母被消去，现在我们只需要计算系统组态的前后变化。上述马尔可夫过程通常称为马尔可夫链。上式显然并不唯一地规定W()，W的显式形式的选择仍然有某些任意性。两种经常使用的选择是： (2.7.10)或 (2.7.11)其中s是一个任意因子，暂时可以取它为1。因此上式可写为： (2.7.12)上式表示能量增加的状态跃迁几率小于1，能量减少的状态跃迁几率为1。因此用上式可以确定马尔可夫链进行的方向，当H0时，状态由，若H0时，则进一步用下述进行判断，e

6、xp(-H) 不允许exp(-H) 允许其中值是一个0与1之间的随机数值。上式说明能量增加较少的状态是允许的，能量增加较多的状态不允许，判断标准是用随机数，物理上原因是体系存在涨落或测不准关系。在马尔可夫链经过N步（N很大）以后，可以认为体系从随机的初始状态出发最后达到了平衡态的附近，对以后的马尔可夫链上的微观状态我们继续按马尔可夫链抽样计算平均值 (2.7.13)相应的偏差计算公式为(2.7.14)N+MNHH马尔可夫初态（随机）上面的马尔可夫过程，对能量的情况，画出图解如下图：图3：能量H的马尔可夫链数学上可以证明，只要N足够大，马尔可夫链出现的状态与初态无关。最后总是到达平衡态的附近。(

7、2.7.12)式子给出的跃迁概率可以使得各状态x0，x1，xn的分布最终为Boltzmann分布。可以直观地论证如下。设我们仅限于处理只发生单个自旋Si反转的试图：假设把自旋Si反转为-Si时会损失能量。由于我们总是想要处于或者靠近模型的基态，所以我们应当以概率为1接受这一变动。因此，在能量变化H=H(xnew)-H(xold)为负的情形下，我们有W(xnew|xold)=1。但是，这样一来我们肯定会陷在一个局部能量极小中。为了避免发生这种情况，我们也得接受使能量增加的变动。但是，我们只允许使能量增加很多的变动很少地发生，因此它们发生的概率应当很低。反之，如果能量变化很小，即原位形同新位形的能

8、量接近，那么我们允许这种变动以相当高的概率发生。这样我们就能向上攀登，挣脱局部能量极小。(2.7.12)式子给出的跃迁概率是一种直观上看来合理、而且也可以精确证明的，它叫做Metropolis函数。但是，这并不是对跃迁概率的唯一可能的选取，还有别的选取法。（二）随机模拟的Monte-Carlo方法在磁性材料模拟计算中的具体应用：这里模拟了单层铁磁和铁磁/反铁磁双层膜系统，一切自旋初始的指向为随机分布，对于薄膜层中选择一个格点调用随机函数进行自旋旋转，计算这个旋转相关的能量改变，同时计算这一旋转的跃迁几率或归一化跃迁几率，产生0，1区间的随机数，再进行进一步的判断，对于满足判定条件的自旋进行旋转

9、。我们建立了单个自旋反转的跃迁概率。一旦对所有的自旋都给过一次反转方向的机会，就完成了一次扫描。一次扫描也叫做每个自旋的一个蒙特卡罗步，缩写为MCS。重要抽样并不是一开始就能运作的，我们利用一个有序位形来启动它，然后源源不断生成新状态。由于两个状态只相差一个自旋转向，因此它们的物理性质有很强的关联，可以经过n个MCS计算一次模拟参数。模拟程序设计具体的模拟是通过对算法的透彻理解，采用Fortran语言编写出相应的程序（详见附录），通过对程序的调试，再根据结果对程序作出相应的修改，最终得出完善的程序。以下是程序提纲。（一）初始化函数调用初始化函数initial，对系统进行初始化，使铁磁层粒子自旋

10、方向均沿着易轴方向，这里设易轴方向沿Z轴正向，此时自旋能量为最大值1；使反铁磁层粒子按照需要取不同的自旋状态，与此同时，在初始函数中可以对反铁磁层的掺杂浓度和方式进行调节。（二）场冷却当系统存在反铁磁层时调用冷却场子函数cool，对系统进行场冷却。1、这里我们对反铁磁层中，没有掺杂的格点进行冷却，并且在平面即X、Y方向采用周期性边界条件，在垂直方向即Z方向采用自由边界条件。2、调用状态函数new，产生系统的新状态。如果新状态下的粒子自旋能量小于原来系统能量，则粒子自旋按要求进行翻转或偏转；如果新状态下的粒子自旋能不小于原状态，那么则用下式进一步判断：exp(-H) 不允许exp(-H) 允许其

11、中值是一个0与1之间的随机数值。上式说明能量增加较少的状态是允许的，能量增加较多的状态不允许，判断标准是用随机数，物理上原因是体系存在涨落或测不准关系。3、由于马尔可夫链要经历一段时间才能使系统状态接近于我们所要的平衡分布，因此我们将一次模拟初始尚未很好地达到平衡分布的那些状态排除在最后分析之外。在本次工作中，一共扫描10000个蒙特卡罗步（MCS），将前面的5000个舍去，然后在剩下的5000步中，每隔20个MCS做一次记录，然后进行统计平均。在一个外加磁场下，求得一个系统的平均磁化强度，从而得出磁化强度随外磁场的变化关系。（三）磁滞回线的测量1、外场为-0.20.2，步长为0.008的条件

12、下测量磁滞回线，同时记录各个点铁磁层的磁化畴结构。2、调用热浴子函数heatbath对系统进行热浴，同时边测量边输出。3、调用能量子函数，利用系统哈密顿量公式计算出系统的能量，求出相应的交换偏置和矫顽场。附：模拟磁矩翻转机制程序 !磁矩翻转机制研究dz=0.01 dy=0.01 jfm=1.0 单晶FM xy格点（三角形材料）! program dilution !单晶单层纯铁磁膜!计算次近邻及次次近邻磁偶极相互作用 parameter(ngx=109,ngy=109) !材料规格 dimension mark(ngx,ngy),spinx(ngx,ngy,2) dimension spiny

13、(ngx,ngy,2),spinz(ngx,ngy,2) dimension zangle(ngx,ngy),yangle(ngx,ngy) dimension ry(ngx,ngy,1),rz(ngx,ngy,1) integer imcs,nmcs,i,x,y,k,l,l1,ngdx,ngdy integer nx,ny,ndegree real kt,b,m,u,v,q,w,dx,dy,dz,bc,Jfm real fi,vxy dimension dsig(25) !下降支取样点标识 dimension usig(25) !上升支取样点标识 integer ksig !取样角度标识（4

14、个角度） j=0 !取样起始位置 do i=1,25 dsig(i)=j usig(i)=j j=j+1 enddo nx=7 !起始格点数x ny=7 !起始格点数y ngdx=103 !末格点数x ngdy=103 !末格点数y len=7 kt=0.1 ! 初始温度 Jfm=1.0 !铁磁的交换耦合系数 bc=0 dx=0 ! 铁磁x方向的单向各向异性常数 dy=0 ! 铁磁y方向的单向各向异性常数 dz=0.05 ! 铁磁z方向的单向各向异性常数 nmcs=2 !蒙特卡诺步 pi=3.1415926 ksig=0 ksig=ksig+1 l=0 l1=0 v=0.00555556 !

15、外场与z方向的角度为1度 ndegree=1 u=v*PI q=cos(u) w=sin(u) call initial(len,mark,spinx,spiny,spinz,zangle, &yangle,ngdx,ngdy,nx,ny,ngx,ngy) do b=0.5,-0.5,-0.04 l=l+1 do imcs=1,nmcs call heatbath(len,kt,b,mark,spinx,spiny,spinz,zangle,yangle,q,w &,dx,dy,dz,Jfm,ngdx,ngdy,nx,ny,ngx,ngy) end do call output(b,spinz

16、,spiny,q,w,k,ngdx,ngdy,nx,ny,ngx,ngy,ndegree) !if(ksig=1.or.ksig=2.or.ksig=3.or.ksig=4)then ! 取样角度为 0 10 60 70 if(.true.)then do i=1,12 if(l=dsig(i)then call keepstatedown(pi,spinx,spiny,spinz,k,l,dsig,i, &ngdx,ngdy,nx,ny,ngx,ngy,ndegree) endif end do end if end do do b=-0.5,0.5,0.04 l1=l1+1 do imcs

17、=1,nmcs call heatbath(len,kt,b,mark,spinx,spiny,spinz,zangle,yangle,q,w &,dx,dy,dz,Jfm,ngdx,ngdy,nx,ny,ngx,ngy) end do call output(b,spinz,spiny,q,w,k,ngdx,ngdy,nx,ny,ngx,ngy,ndegree) !if(ksig=1.or.ksig=2.or.ksig=3.or.ksig=4)then if(.true.)then do i=1,12 if(l1=usig(i)then call keepstateup(pi,spinx,s

18、piny,spinz,k,l1,usig,i, &ngdx,ngdy,nx,ny,ngx,ngy,ndegree) endif end do end if end do end program subroutine initial(len,mark,spinx,spiny,spinz,zangle,yangle, &ngdx,ngdy,nx,ny,ngx,ngy) dimension mark(ngx,ngy),spinx(ngx,ngy,2) dimension spiny(ngx,ngy,2),spinz(ngx,ngy,2) dimension eangle(ngx,ngy),zangl

19、e(ngx,ngy),yangle(ngx,ngy) real newx,newy,newz integer x,y,k pi=3.1415926 do y=1,ngx do x=1,ngy mark(x,y)=1 spinx(x,y,1)=0 spiny(x,y,1)=0 spinz(x,y,1)=0 enddo enddo do y=ny,ngdy do x=nx,ngdx if(4*y+7*x-721)0.and.(4*y-7*x+49)0)then !三角形 mark(x,y)=0 spinx(x,y,1)=0 spiny(x,y,1)=0 spinz(x,y,1)=1 else ma

20、rk(x,y)=1 spinx(x,y,1)=0 spiny(x,y,1)=0 spinz(x,y,1)=0 endif end do end do end subroutine heatbath(len,kt,b,mark,spinx,spiny,spinz,zangle, &yangle,q,w,dx,dy,dz,Jfm,ngdx,ngdy,nx,ny,ngx,ngy) dimension mark(ngx,ngy),spinx(ngx,ngy,2) dimension spiny(ngx,ngy,2),spinz(ngx,ngy,2) dimension zangle(ngx,ngy),

21、yangle(ngx,ngy) integer x,y,z real ee,newx,newy,newz,kt,q,w,Jfm z=1 do x=nx,ngdx do y=ny,ngdy if(mark(x,y)/=1)then call e(len,x,y,z,b,ee,newx,newy,newz,mark,spinx,spiny,spinz, &zangle,yangle,q,w,dx,dy,dz,Jfm,ngdx,ngdy,nx,ny,ngx,ngy) call RANECU(RVEC,LEN) if(ee0.or.RVEC=1/(1+exp(ee/kt)then ! 判断能量改变是否

22、被允许 spinx(x,y,z)=newx spiny(x,y,z)=newy spinz(x,y,z)=newz end if endif end do end do end subroutine e(len,x,y,z,b,ee,newx,newy,newz,mark,spinx,spiny, &spinz,zangle,yangle,q,w,dx,dy,dz,Jfm,ngdx,ngdy,nx,ny,ngx,ngy) dimension mark(ngx,ngy),spinx(ngx,ngy,2) dimension spiny(ngx,ngy,2),spinz(ngx,ngy,2) di

23、mension zangle(ngx,ngy),yangle(ngx,ngy) real leftx,rightx,behindx,frontx,lefty,righty,behindy,fronty,dz real leftz,rightz,behindz,frontz,newx,newy,newz,Jint,Jafm,Jfm,kz real dy real EDD1,EDD2 real sumx,sumy,sumz,e1x,e1k,e1z,e2x1,e2x2,e2k,e2z,e1,e2,ee,q,w integer x,y,z Jint=0 if(x=nx)then leftx=0 ! 畴

24、壁间无相互作用 lefty=0 leftz=0 else ! 计算畴结构中当前格点与临近格点相互作用中临近格点的等效 !自旋 leftx=spinx(x-1,y,z) lefty=spiny(x-1,y,z) leftz=spinz(x-1,y,z) end if if(x=ngdx)then rightx=0 righty=0 rightz=0 else rightx=spinx(x+1,y,z) righty=spiny(x+1,y,z) rightz=spinz(x+1,y,z) end if if(y=ny)then behindx=0 behindy=0 behindz=0 else

25、 behindx=spinx(x,y-1,z) behindy=spiny(x,y-1,z) behindz=spinz(x,y-1,z) end if if(y=ngdy)then frontx=0 fronty=0 frontz=0 else frontx=spinx(x,y+1,z) fronty=spiny(x,y+1,z) frontz=spinz(x,y+1,z) end if call new(len,newx,newy,newz) ! 产生新的自旋 sumx=leftx+rightx+behindx+frontx ! 临近格点的等效自旋x方向分量 sumy=lefty+righ

26、ty+behindy+fronty sumz=leftz+rightz+behindz+frontz call DD(x,y,z,b,ee,newx,newy,newz,mark,spinx,spiny,spinz,zangle &,yangle,q,w,dx,dy,dz,Jfm,ngdx,ngdy,EDD1,EDD2,nx,ny,ngx,ngy) if(z=1)then e1x=-Jint*(newx*spinx(x,y,2)+newy*spiny(x,y,2)+newz*spinz(x,y,2) ! 铁磁层自旋改变后交换耦合能 e1k=-dz*newz*2-dy*newy*2 ! 铁磁层自

27、旋改变后各向异性能 e1z=-b*newz*q-b*newy*w-Jfm*(newx*sumx+newy*sumy+newz*sumz) ! 自旋改变后塞曼能和铁磁交换耦合能 e1=e1x+e1k+e1z+EDD1 ! 铁磁层末能 e2x=-Jint*(spinx(x,y,z)*spinx(x,y,2)+spiny(x,y,z) &*spiny(x,y,2)+spinz(x,y,z)*spinz(x,y,2) ! 自旋改变前界面耦合能 e2k=-dz*(spinz(x,y,z)*2-dy*(spiny(x,y,z)*2 ! 铁磁层自旋改变前各向异性能 e2z=-b*spinz(x,y,z)*q

28、-b*spiny(x,y,z)*w-Jfm*(spinx(x,y,z) &*sumx+spiny(x,y,z)*sumy+spinz(x,y,z)*sumz) ! 铁磁层塞曼能和交换耦合能 e2=e2x+e2k+e2z+EDD2 ! 铁磁层初能 ee=e1-e2 ! print*,e1 !print*,e2 else e1x=-Jint*(newx*spinx(x,y,1)+newy*spiny(x,y,1)+newz*spinz(x,y,1) ! 铁磁层自旋改变后交换耦合能 e1k=-dz*newz*2-dy*newy*2 ! 铁磁层自旋改变后各向异性能 e1z=-b*newz*q-b*ne

29、wy*w-Jfm*(newx*sumx+newy*sumy+newz*sumz) ! 自旋改变后塞曼能和铁磁交换耦合能 e1=e1x+e1k+e1z+EDD1 ! 铁磁层末能 e2x=-Jint*(spinx(x,y,z)*spinx(x,y,1)+spiny(x,y,z) &*spiny(x,y,1)+spinz(x,y,z)*spinz(x,y,1) ! 自旋改变前界面耦合能 e2k=-dz*(spinz(x,y,z)*2-dy*(spiny(x,y,z)*2 ! 铁磁层自旋改变前各向异性能 e2z=-b*spinz(x,y,z)*q-b*spiny(x,y,z)*w-Jfm*(spinx

30、(x,y,z) &*sumx+spiny(x,y,z)*sumy+spinz(x,y,z)*sumz) ! 铁磁层塞曼能和交换耦合能 e2=e2x+e2k+e2z+EDD2 ! 铁磁层初能 ee=e1-e2 end if end subroutine DD(x,y,z,b,ee,newx,newy,newz,mark,spinx,spiny,spinz, &zangle,yangle,q,w,dx,dy,dz,Jfm,ngdx,ngdy,EDD1,EDD2,nx,ny,ngx,ngy) integer x1,y1,z1,x,y,z dimension mark(ngx,ngy),spinx(ngx,ngy,2) dimension spiny(ngx,ngy,2),spinz(ngx,ngy,2) dimension zangle(ngx,ngy),yangle(ngx,ngy) real ee