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

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

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

3、型，各個(gè)自旋傾向于反向排列。2.模擬二維Ising模型的步驟：為了方便，令為1。(1) 選擇任意一個(gè)初始位形Xx1,x2, xN；(2) 按1/N等概率的選取一個(gè)格點(diǎn)i，將其自旋反向，得到一個(gè) 新的位形Xx1,x2, xN；(3) 利用公式 ，計(jì)算能量差E=E(X)-E(X) ，若E<0，則改變有效，位形改變X X；(4) 如果E>0，則再產(chǎn)生一個(gè)0,1之間的隨機(jī)數(shù)，如果 ；則位形改變有效，否則位形不變；(5) 返回步驟(2)，進(jìn)行下一次迭代。（）固定外加磁場(chǎng)，磁化強(qiáng)度隨溫度的變化源程序#include <iostream> #include <string.h&

4、gt;#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 randx,randnum;unsigned long randxi=1;const unsigned long randa=16807;unsigne

5、d 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; fp8=fopen("固定外加磁場(chǎng)，隨溫度.txt","w+"); const int le

6、nx=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,ipy=0; double H=0.0; for (i=0;i<lenx;i+) imi=i-1; ipi=i+1

7、; im0=lenx; iplenx-1=0; /設(shè)置邊界條件 mag0=lenx*leny; magmax=abs(mag0);/磁矩最大絕對(duì)值 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;iy<leny;iy+) spinixiy=1;/初始狀態(tài)，設(shè)置所有初始磁矩為1

8、 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;/ 對(duì)格點(diǎn)（ix，iy）周圍最近鄰的磁矩求和 den=2.0*enj*spinkxky*dspin+2*H*spinkxky; /能量變化，由于磁矩從-1到1或1到-1，所以乘系數(shù)2 trapro=ex

9、p(-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"<<'n' goto labelout; la

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

11、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=2147483637;/ random number creation *void random_number()

12、 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("固定溫度，隨外加磁場(chǎng).txt","w+"); const int lenx=100,leny=lenx; int spinlenxleny,imlenx,iplenx; double mag0=0.0,mag=

13、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 for (i=0;i<lenx;i+) imi=i-

14、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;/磁場(chǎng)從-10增加到10 else h=h-0.02; /磁場(chǎng)從10減少到-10 for (mcsi=0;mcsi<m

15、csmax;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;/ 對(duì)格點(diǎn)（ix，iy）周圍最近鄰的磁矩求和 den=2.0*enj*spinkxky*dspin+2*h*spinkxky;/能量變化，由于磁矩從-1到1或1到-1，所以乘系數(shù)2 trapro=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,"%