二維Ising模型的程序設(shè)計(jì)

1、二維Ising模型的程序設(shè)計(jì)一、課題名稱：二維Ising模型的程序設(shè)計(jì)二、班級和姓名：*三、主要內(nèi)容：1.研究的內(nèi)容和算法：Ising模型最初由Lenz提出和用來作為鐵磁性的一個(gè)模型。后來 成為他的研究生Ising的博士論文的題目。1925年，Ising給出了一維 情況下的解，該解顯示，在一維情況下，Ising模型沒有相變解。1944 年，Onsager得到了二維Ising模型的準(zhǔn)確解，二維時(shí)就有了相變。 對于三維，至今還沒有嚴(yán)格解，需依靠數(shù)值計(jì)算得到。物質(zhì)在外磁場H中的磁場強(qiáng)度乂為“ H(1)抗磁體，x<0,數(shù)值很小且是常數(shù)，不隨溫度變化；(2)順磁體，>>0,數(shù)值很小且隨

2、溫度反比或與溫度無關(guān)；(3)鐵磁體，在一定相變溫度 Tc (Curie溫度)之下，M不隨 H作線性變化，具有磁滯回線是磁體物質(zhì)的在磁場中行為的基本特 性，磁化率與外磁場有關(guān)。在 Tc之上時(shí)，鐵磁性消失，轉(zhuǎn)變?yōu)轫槾?性。(4)反鐵磁體，溫度在Tc之上時(shí)是順磁體，之下時(shí)x隨溫度下 降而降低。對于二維Ising模型，令：G=Ld為一個(gè)d維、共有N個(gè)格點(diǎn)的 體系，在每個(gè)格點(diǎn)i上有一個(gè)自旋，可以朝上或朝下的方向。用自旋變量o-i表示,1,自旋朝上1,自旋朝下在外磁場H中，體系的哈密頓量為：2i1 i，ji 1 ,其中J為交換關(guān)聯(lián)系數(shù)，pB表示單個(gè)自旋的磁矩，<i,j>表示只對格 點(diǎn)i周圍最鄰

3、近的給點(diǎn)j求和。J為正時(shí)為鐵磁體的模型，各個(gè)自旋傾向于同向排列；J為負(fù)時(shí)為反磁體的模型，各個(gè)自旋傾向于反向排列。2模擬二維Ising模型的步驟：為了方便，令”,訪1。(1)選擇任意一個(gè)初始位形 Xx1,x2xN;(2)按1/N等概率的選取一個(gè)格點(diǎn)i,將其自旋反向，得到一個(gè)新的位形 X'x1,x2xN;NNH Ei j B H i 利用公式2 i 1 i,j ji 1,計(jì)算能量差E=E(X)-E(X),若E<0,則改變有效，位形改變 X- X' 如果 E>0 ,則再產(chǎn)生一個(gè)0,1之間的隨機(jī)數(shù) J 如 果 e E ；則位形改變有效，否則位形不變；(5)返回步驟(2),進(jìn)

4、行下一次迭代。(I)固定外加磁場，磁化強(qiáng)度隨溫度的變化源程序#include <iostream>#include <string.h>#include <math.h>#include <iomanip>#include <stdlib.h>#include <fstream>#include <cstdio>#include <iomanip>using namespace std;const double pi=3.1415926;const double em=1.0e-12;double

5、randx,randnum;unsigned long randxi=1;const unsigned long randa=16807;*unsigned long randm=2147483637;/ random number creation void random_number()if ( randxi = randm )randxi=randm-1;randxi=(randa*randxi)%randm;randx=randxi;randnum=randx/randm;return;/end random_number*/ main*int main(void)FILE *fp8;

6、fp8=fopen("固定外加磁場，隨溫度.txt","w+");const int lenx=100,leny=lenx;int spinlenxleny,imlenx,iplenx;double mag0=0.0,mag=0.0;double enj=0.5,magmax=0.0;double temp=0.0,den=0.0,trapro=0.0;int l=0,mcsi=0,mcsmax=0;int i=0,j=0,ix=0,iy=0,iz=0,kx=0,ky=0,kz=0,t=0;int dspin=0,imx=0,ipx=0,imy=0,ip

7、y=0;double H=0.0;for (i=0;i<lenx;i+)imi=i-1;ipi=i+1;im0=lenx;iplenx-1=0;/設(shè)置邊界條件mag0=lenx*leny;magmax=abs(mag0);/磁矩最大絕對值cout<<setw(15)<<"H= "<< 'n'cin>>H;mcsmax=pow(10,5);for (i=1;i<=100;i=i+1) temp=temp+0.1;mag=mag0;for (ix=0;ix<lenx;ix+)for (iy=0;i

8、y<leny;iy+)spinixiy=1;初始狀態(tài)，設(shè)置所有初始磁矩為for (mcsi=0;mcsi<mcsmax;mcsi+)random_number();kx=int(lenx*randnum);random_number();ky=int(leny*randnum);imx=imkx;ipx=ipkx;imy=imky;ipy=ipky；dspin=spinimxky+spinipxky+spinkximy+spinkxipy;/對格點(diǎn)(ix, iy)周圍最近鄰的磁矩求和den=2.0*enj*spinkxky*dspin+2*H*spinkxky;/ 能量變化，由于磁

9、矩從-1到1或1至U-1 ,所以乘系數(shù)2trapro=exp(-den/temp);random_number();if (den <= 0 | randnum <= trapro) mag=mag-2.0*spinkxky;/ 磁矩的變化spinkxky=-spinkxky; / 自旋反轉(zhuǎn) else goto label1;label1: continue; fprintf(fp8,"%15.6f %15.6f n", temp,mag/magmax);cout<<setw(10)<<"OK"<<'

10、;n'goto labelout;labelout: return 0;/ end main*計(jì)算結(jié)果及具體分析討論取關(guān)聯(lián)系數(shù)為0.5,外加磁場分別取0.5, 1.0, 1.5,得到磁化強(qiáng)度隨溫度的變化關(guān)系如下圖:OJOjGLI L6，i 2.53.13石 4,14.6519.6由圖可得，在不同的外加磁場下，隨著溫度的增加，磁化強(qiáng)度降低，這應(yīng)該是溫度增加導(dǎo)致的電子的熱運(yùn)動增加，電子運(yùn)動變得更加無序，所以磁矩變小。(II)固定溫度，磁化強(qiáng)度隨外加磁場的變化源程序#include <iostream>#include <string.h>#include <m

11、ath.h>#include <iomanip>#include <stdlib.h>#include <fstream>#include <cstdio>#include <iomanip>using namespace std;const double pi=3.1415926;const double em=1.0e-12;double randx,randnum;unsigned long randxi=1;const unsigned long randa=16807;*unsigned long randm=2147

12、483637;/ random number creation void random_number()if ( randxi = randm )randxi=randm-1;randxi=(randa*randxi)%randm;randx=randxi;randnum=randx/randm;return;/end random_number*/ main*int main(void)FILE *fp8;fp8=fopen("固定溫度，隨外加磁場.txt","w+");const int lenx=100,leny=lenx;int spinlenx

13、leny,imlenx,iplenx;double mag0=0.0,mag=0.0;double enj=0.9,magmax=0.0;double temp=0.0,den=0.0,trapro=0.0;int l=0,mcsi=0,mcsmax=0;int i=0,j=0,ix=0,iy=0,iz=0,kx=0,ky=0,kz=0,t=0;int dspin=0,imx=0,ipx=0,imy=0,ipy=0;double h=-10;for (ix=0;ix<lenx;ix+)for (iy=0;iy<leny;iy+)spinixiy=1;/初始狀態(tài)，設(shè)置所有初始磁矩為1

14、)for (i=0;i<lenx;i+)(imi=i-1;ipi=i+1；)im0=lenx;iplenx-1=0;mag0=lenx*leny;magmax=abs(mag0);mag=mag0;cout<<setw(15)<<"Temp= "<< 'n'cin>>temp;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減少到-10f

15、or (mcsi=0;mcsi<mcsmax;mcsi+)(random_number();kx=int(lenx*randnum); random_number();ky=int(leny*randnum);imx=imkx;ipx=ipkx;imy=imky;對格點(diǎn)（ix,ipy=ipky;dspin=spinimxky+spinipxky+spinkximy+spinkxipy;/ iy)周圍最近鄰的磁矩求和den=2.0*enj*spinkxky*dspin+2*h*spinkxky;/ 能量變化，由于磁矩從-1到1或1到-1 ,所以乘系數(shù)2trapro=exp(-den/temp);random_number();if (den <= 0 | randnum <= trapro)mag=mag-2.0*spinkxky;/磁矩的變化spinkxky=-spinkxky;/ 自旋反轉(zhuǎn)) else goto label1;) labell: continue;fprintf(fp8,"%15.6f %