[原创]用Monte-Carlo方法研究二维Ising模型的相变问题，求出... - ChinaNet

/**************************************************************************************
*程序用途：用Monte-Carlo方法研究二维Ising模型的相变问题，求出临界温度T，比热C及磁化率χ
*编译环境：Visual C++ 6.0
*作者：lxg
*单位：
*完成时间：2005年12月28日
**************************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define L 150                      /*L×L个格点                         */
#define preN 100000                /*舍弃前preN个E，M值                */
#define lateN 5100000                /*总共计算lateN个E，M值             */
#define deltaN (lateN-preN)          /*Monte-Carlo方法抽样E，M个数       */
double recordEnergy[deltaN];       /*存储抽样能量值                   */
int recordMagnet[deltaN];          /*存储抽样磁矩量                   */
int lattice[L][L];                   /*格点，存储格点状态                */
int temp[L][L];                      /*保留初始化格点状态                */
double energy;                      /*记录初试化及试探后系统能量，磁矩 */
int magnet;                         /*并赋给recordEnergy和recordMagnet */
double T;                            /*温度                               */
struct CX                            /*比热C和磁化率χ                   */
{
double C;
double X;
};

/****************************************************************************
*函数原型
*****************************************************************************/
void InitialLattice();
void SaveInitial();
void GetInitial();
void CaculateFirstMagnet();
void CaculateFirstEnergy();
void GenerateNext();
struct CX Caculate_C_and_X();

/****************************************************************************
*主函数
*****************************************************************************/
int main()
{
int i;
struct CX cx[200];
FILE *ft,*fx,*fc;
ft=fopen("T.txt","w");
fx=fopen("X.txt","w");
fc=fopen("C.txt","w");
srand(time(NULL));

T = 1.0;
for(i=0;T<4.0;i++)
{
  InitialLattice();
  cx=Caculate_C_and_X();
  fprintf(ft,"%f\n",T);
  fprintf(fx,"%f\n",cx.X);
  fprintf(fc,"%f\n",cx.C);
  printf("T=%f ",T);
  printf("  X=%f ",cx.X);
  printf("  C=%f\n",cx.C);
      T = T+0.05;
}

fclose(ft);
fclose(fx);
fclose(fc);
return 0;
}

/***************************************************************************
*函数实现
****************************************************************************/
void InitialLattice()                            /*初始化格点状态，可以  */
{                                                 /*通过热初始化和冷初始化*/
int i,j;
for(i=0;i<L;i++)
{
      for(j=0;j<L;j++)
  {
if((rand()%10)>5) lattice[j]=1;
else             lattice[j]=1;    /*1冷初始化，-1热初试化*/
      }
}
}

void SaveInitial()                               /*保存第一次初始化状态*/
{                                                 /*若为冷初试化，可以不*/
int k,m;                                     /*调用该函数       */
for(k=0;k<L;k++)
{
  for(m=0;m<L;m++)
temp[k][m]=lattice[k][m];
}
}

void GetInitial()                                  /*提取第一次初始化状态*/
{                                                 /*若为冷初试化，可以不*/
int k,m;                                     /*调用该函数       */
for(k=0;k<L;k++)
{
  for(m=0;m<L;m++)
lattice[k][m]=temp[k][m];

}
}

void CaculateFirstEnergy()                         /*计算初试化后系统能量*/
{

double dEnergy = 0.0;
int up,down,right,left,px,py;
for(px=0;px<L;px++)
{
  for(py=0;py<L;py++)
  {
up = (L+px-1)%L;
down = (px+1)%L;
right = (L+py-1)%L;
left = (py+1)%L;
dEnergy += 0.5*lattice[px][py]*(lattice[up][py]+lattice[down][py]+lattice[px][right]+lattice[px][left]);
  }
}

energy = dEnergy;
}

void GenerateNext()                                  /*试探，产生下一状态*/
{
int px,py;
int dE,sum;
int up,down,right,left;
px = rand()%L;
py = rand()%L;

up = (L+px-1)%L;
down = (px+1)%L;
right = (L+py-1)%L;
left = (py+1)%L;

sum = (lattice[px][py])*(lattice[up][py]+lattice[down][py]+lattice[px][right]+lattice[px][left]);
switch(sum)
{
case 4: dE= 8; break;
case 2: dE= 4; break;
case 0: dE= 0; break;
case -2: dE= -4; break;
case -4: dE= -8; break;
default: printf("---Erro!---\n"); break;
}

if( dE>0 )                                           /*判断是否接受试探*/
{
      if((rand()*1.0/RAND_MAX) <= exp((-1.0)*dE/T))
  {
lattice[px][py] = -lattice[px][py];
magnet += 2*lattice[px][py];                   /*更新磁矩大小*/
energy += dE;                                  /*更新能量大小*/
  }

}
else
{
  lattice[px][py] = -lattice[px][py];
  magnet += 2*lattice[px][py];                      /*更新磁矩大小*/
  energy += dE;                                     /*更新能量大小*/
}

}

struct CX Caculate_C_and_X()                               /*计算比热C及磁化率χ*/
{
struct CX cx;
int i;
double AvgEnergy=0.0;
double AvgSquareEnergy=0.0;
double AvgMagnet=0.0;
double AvgSquareMagnet=0.0;

CaculateFirstMagnet();                                  /*初始化磁矩*/
CaculateFirstEnergy();                                  /*初始化能量*/

for(i=0; i<lateN; i++)
{
      GenerateNext();                                     /*试探       */
      if(i>=preN)                                        /*舍弃前preN个*/
      {
         recordEnergy[i-preN] = energy;                   /*记录能量值  */
         recordMagnet[i-preN] = magnet;                   /*记录磁矩值  */
  }
}

for(i=0;i<deltaN;i++)                                  /*计算物理量的平均值*/
{
      AvgEnergy += (double)recordEnergy/deltaN;
      AvgSquareEnergy += (double)recordEnergy*recordEnergy/deltaN;
      AvgMagnet += (double)recordMagnet/deltaN;
      AvgSquareMagnet += (double)recordMagnet*recordMagnet/deltaN;
}

cx.C = (AvgSquareEnergy - AvgEnergy*AvgEnergy)/T/T;       /*得出结果，并返回*/
cx.X = (AvgSquareMagnet - AvgMagnet*AvgMagnet)/T;

return cx;
}