用Monte-Carlo方法研究二維Ising模型的相變問題-Read

1、 原創(chuàng) 用 Monte-Carlo 方法研究二維 Ising 模型的相變問題， 求出臨界溫度 T，比熱 C 及磁化率 圖片：圖片：/*程序用途：用 Monte-Carlo 方法研究二維 Ising 模型的相變問題，求出臨界溫度T，比熱 C及磁化率*編譯環(huán)境： Visual C+ 6.0*作 者：lxg* 單 位：*完成時(shí)間： 2005 年 12 月 28 日*/#include#include#include#include#define L 150#define preN 100000#define lateN 5100000#define deltaN (lateN-preN)/*L L

2、個(gè)格點(diǎn)*/* 舍棄前 preN 個(gè) E， M 值*/*總共計(jì)算 lateN 個(gè) E， M 值*/*Monte-Carlo 方法抽樣 E ，M 個(gè)數(shù) */double recordEnergydeltaN; int recordMagnetdeltaN;int latticeLL;int tempLL;double energy;int magnet;double T;struct CXdouble C;double X;/* 存儲(chǔ)抽樣能量值/* 存儲(chǔ)抽樣磁矩量/*格點(diǎn)，存儲(chǔ)格點(diǎn)狀態(tài)/*保留初始化格點(diǎn)狀態(tài)/*記錄初試化及試探后系統(tǒng)能量，磁矩*/*/*/* 并賦給 recordEnergy 和 r

3、ecordMagnet /* 溫度*/* 比熱 C 和磁化率 */*/*/*/* *函數(shù)原型*/ struct CX Caculate_C_and_X();void voidvoidvoid void voidInitialLattice();SaveInitial();GetInitial();CaculateFirstMagnet();CaculateFirstEnergy();GenerateNext();/* *主函數(shù)*/ int main()int i;struct CX cx200; FILE *ft,*fx,*fc; ft=fopen(T.txt,w); fx=fopen(X.t

4、xt,w); fc=fopen(C.txt,w);srand(time(NULL);T = 1.0;for(i=0;T4.0;i+)InitialLattice();cx =Caculate_C_and_X(); fprintf(ft,%fn,T); fprintf(fx,%fn,cx.X); fprintf(fc,%fn,cx.C);printf(T=%f ,T); printf( X=%f ,cx.X); printf( C=%fn,cx.C);T = T+0.05;fclose(ft); fclose(fx); fclose(fc); return 0; /* 保存第一次初始化狀態(tài) */

5、 /* 若為冷初試化，可以不 */ /*調(diào)用該函數(shù)*/* 提取第一次初始化狀態(tài) */ /* 若為冷初試化，可以不 */ /*調(diào)用該函數(shù)*/*函數(shù)實(shí)現(xiàn) */ void InitialLattice()/*初始化格點(diǎn)狀態(tài)，可以 */* 通過熱初始化和冷初始化 */int i,j;for(i=0;iL;i+)for(j=0;j5) latticej=1;else latticej=1;/*1 冷初始化， -1 熱初試化 */void SaveInitial()int k,m;for(k=0;kL;k+)for(m=0;mL;m+) tempkm=latticekm;void GetInitial()

6、int k,m;for(k=0;kL;k+)for(m=0;mL;m+)latticekm=tempkm;void CaculateFirstMagnet()/* 計(jì)算初試化后系統(tǒng)磁矩 */int i,j;int iMagnet=0;for(i=0;iL;i+) for(j=0;jL;j+) iMagnet += latticej;magnet = iMagnet;void CaculateFirstEnergy()/* 計(jì)算初試化后系統(tǒng)能量 */double dEnergy = 0.0;int up,down,right,left,px,py;for(px=0;pxL;px+)for(py=

7、0;py0 )/*判斷是否接受試探 */if(rand()*1.0/RAND_MAX) = exp(-1.0)*dE/T)latticepxpy = -latticepxpy; magnet += 2*latticepxpy; energy += dE;/* 更新磁矩大小 */* 更新能量大小 */elselatticepxpy = -latticepxpy; magnet += 2*latticepxpy; energy += dE;/*更新磁矩大小 */ /* 更新能量大小 */struct CX Caculate_C_and_X()struct CX cx;int i;double Av

8、gEnergy=0.0;double AvgSquareEnergy=0.0;double AvgMagnet=0.0;double AvgSquareMagnet=0.0;/*計(jì)算比熱 C及磁化率 */*初始化磁矩 */*初始化能量 */CaculateFirstMagnet();CaculateFirstEnergy();for(i=0; i=preN)recordEnergyi-preN/*舍棄前preN 個(gè) */= energy;/* 記錄能量值 */recordMagneti-preN= magnet;/* 記錄磁矩值 */for(i=0;ideltaN;i+)/* 計(jì)算物理量的平均值 */AvgEnergy += (double)recordEnergy/deltaN;AvgSquareEnergy += (double)recordEnergy*recordEnergy/deltaN;AvgMagnet += (double)recordMagnet/deltaN;AvgSquareMagnet += (double)recordMagnet*recordMagnet/deltaN;cx.