二维Ising模型的程序设计

二维Ising模型的程序设计一、课题名称：二维Ising模型的程序设计二、班级和姓名：***三、主要容：1.研究的容和算法：Ising模型最初由Lenz提出和用来作为铁磁性的一个模型。后来成为他的研究生Ising的博士论文的题目。1925年，Ising给出了一维情况下的解，该解显示，在一维情况下，Ising模型没有相变解。1944年，Onsager得到了二维Ising模型的准确解，二维时就有了相变。对于三维，至今还没有严格解，需依靠数值计算得到。物质在外磁场H中的磁场强度M为HM（1）抗磁体，χ0，数值很小且是常数，不随温度变化；（2）顺磁体，χ0，数值很小且随温度反比或与温度无关；（3）铁磁体，在一定相变温度Tc（Curie温度）之下，M不随H作线性变化，具有磁滞回线是磁体物质的在磁场中行为的基本特性，磁化率与外磁场有关。在Tc之上时，铁磁性消失，转变为顺磁性。（4）反铁磁体，温度在Tc之上时是顺磁体，之下时χ随温度下降而降低。对于二维Ising模型，令：G=Ld为一个d维、共有N个格点的体系，在每个格点i上有一个自旋，可以朝上或朝下的方向。用自旋变量σi表示，自旋朝下自旋朝上,1,1i。在外磁场H中，体系的哈密顿量为：NiiBjijNiiEHJH1,12，其中J为交换关联系数，μB表示单个自旋的磁矩，i,j表示只对格点i周围最邻近的给点j求和。J为正时为铁磁体的模型，各个自旋倾向于同向排列；J为负时为反磁体的模型，各个自旋倾向于反向排列。2.模拟二维Ising模型的步骤：为了方便，令BBk,为1。(1)选择任意一个初始位形X{x1,x2,…xN}；(2)按1/N等概率的选取一个格点i，将其自旋反向，得到一个新的位形X’{x1,x2,…xN}；(3)利用公式，计算能量差△E=E(X’)-E(X)，若△E0，则改变有效，位形改变X→X’；(4)如果△E0，则再产生一个[0,1]之间的随机数ξ，如果；则位形改变有效，否则位形不变；(5)返回步骤(2)，进行下一次迭代。（Ⅰ）固定外加磁场，磁化强度随温度的变化①源程序#includeiostream#includestring.h#includemath.h#includeiomanip#includestdlib.hNiiBjijNiiEHJH1,12Ee#includefstream#includecstdio#includeiomanipusingnamespacestd;constdoublepi=3.1415926;constdoubleem=1.0e-12;doublerandx,randnum;unsignedlongrandxi=1;constunsignedlongranda=16807;unsignedlongrandm=2147483637;//randomnumbercreation***********************voidrandom_number(){if(randxi==randm){randxi=randm-1;}randxi=(randa*randxi)%randm;randx=randxi;randnum=randx/randm;return;}//endrandom_number*****************************//main*****************************************intmain(void){FILE*fp8;fp8=fopen(固定外加磁场，随温度.txt,w+);constintlenx=100,leny=lenx;intspin[lenx][leny],im[lenx],ip[lenx];doublemag0=0.0,mag=0.0;doubleenj=0.5,magmax=0.0;doubletemp=0.0,den=0.0,trapro=0.0;intl=0,mcsi=0,mcsmax=0;inti=0,j=0,ix=0,iy=0,iz=0,kx=0,ky=0,kz=0,t=0;intdspin=0,imx=0,ipx=0,imy=0,ipy=0;doubleH=0.0;for(i=0;ilenx;i++){im[i]=i-1;ip[i]=i+1;}im[0]=lenx;ip[lenx-1]=0;//设置边界条件mag0=lenx*leny;magmax=abs(mag0);//磁矩最大绝对值coutsetw(15)H='\n';cinH;mcsmax=pow(10,5);for(i=1;i=100;i=i+1){temp=temp+0.1;mag=mag0;for(ix=0;ixlenx;ix++){for(iy=0;iyleny;iy++){spin[ix][iy]=1;//初始状态，设置所有初始磁矩为1}}for(mcsi=0;mcsimcsmax;mcsi++){random_number();kx=int(lenx*randnum);random_number();ky=int(leny*randnum);imx=im[kx];ipx=ip[kx];imy=im[ky];ipy=ip[ky];dspin=spin[imx][ky]+spin[ipx][ky]+spin[kx][imy]+spin[kx][ipy];//对格点（ix，iy）周围最近邻的磁矩求和den=2.0*enj*spin[kx][ky]*dspin+2*H*spin[kx][ky];//能量变化，由于磁矩从-1到1或1到-1，所以乘系数2trapro=exp(-den/temp);random_number();if(den=0||randnum=trapro){mag=mag-2.0*spin[kx][ky];//磁矩的变化spin[kx][ky]=-spin[kx][ky];//自旋反转}else{gotolabel1;}label1:continue;}fprintf(fp8,%15.6f%15.6f\n,\temp,mag/magmax);}coutsetw(10)OK'\n';gotolabelout;labelout:return0;}//endmain*************************************②计算结果及具体分析讨论取关联系数为0.5，外加磁场分别取0.5，1.0，1.5，得到磁化强度随温度的变化关系如下图：由图可得，在不同的外加磁场下，随着温度的增加，磁化强度降低，这应该是温度增加导致的电子的热运动增加，电子运动变得更加无序，所以磁矩变小。（Ⅱ）固定温度，磁化强度随外加磁场的变化①源程序#includeiostream#includestring.h#includemath.h#includeiomanip#includestdlib.h#includefstream#includecstdio#includeiomanipusingnamespacestd;constdoublepi=3.1415926;constdoubleem=1.0e-12;doublerandx,randnum;外加磁场unsignedlongrandxi=1;constunsignedlongranda=16807;unsignedlongrandm=2147483637;//randomnumbercreation***********************voidrandom_number(){if(randxi==randm){randxi=randm-1;}randxi=(randa*randxi)%randm;randx=randxi;randnum=randx/randm;return;}//endrandom_number*****************************//main*****************************************intmain(void){FILE*fp8;fp8=fopen(固定温度，随外加磁场.txt,w+);constintlenx=100,leny=lenx;intspin[lenx][leny],im[lenx],ip[lenx];doublemag0=0.0,mag=0.0;doubleenj=0.9,magmax=0.0;doubletemp=0.0,den=0.0,trapro=0.0;intl=0,mcsi=0,mcsmax=0;inti=0,j=0,ix=0,iy=0,iz=0,kx=0,ky=0,kz=0,t=0;intdspin=0,imx=0,ipx=0,imy=0,ipy=0;doubleh=-10;for(ix=0;ixlenx;ix++){for(iy=0;iyleny;iy++){spin[ix][iy]=1;//初始状态，设置所有初始磁矩为1}}for(i=0;ilenx;i++){im[i]=i-1;ip[i]=i+1;}im[0]=lenx;ip[lenx-1]=0;mag0=lenx*leny;magmax=abs(mag0);mag=mag0;coutsetw(15)Temp='\n';cintemp;mcsmax=pow(10,5);for(i=1;i=2000;i=i+1){t=t+1;if(t=1000){h=h+0.02;//磁场从-10增加到10}else{h=h-0.02;}//磁场从10减少到-10for(mcsi=0;mcsimcsmax;mcsi++){random_number();kx=int(lenx*randnum);random_number();ky=int(leny*randnum);imx=im[kx];ipx=ip[kx];imy=im[ky];ipy=ip[ky];dspin=spin[imx][ky]+spin[ipx][ky]+spin[kx][imy]+spin[kx][ipy];//对格点（ix，iy）周围最近邻的磁矩求和den=2.0*enj*spin[kx][ky]*dspin+2*h*spin[kx][ky];//能量变化，由于磁矩从-1到1或1到-1，所以乘系数2trapro=exp(-den/temp);random_number();if(den=0||randnum=trapro){mag=mag-2.0*spin[kx][ky];//磁矩的变化spin[kx][ky]=-spin[kx][ky];//自旋反转}else{gotolabel1;}label1:continue;}fprintf(fp8,%15.6f%15.6f\n,\h,mag/magmax);}coutsetw(10)OK'\n';gotolabelout;labelout:return0;}//endmain*************************************②计算结果及具体分析讨论取关联系数为0.9，画出温度分别为0.25，0.5时的磁化强度随外加磁场的变化关系，如下图：温度越低，磁滞回线的面积越大，低温时磁矩更大，所以矫顽力更大，导致磁滞回线面积更大。