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

DoubleX

/**************************************************************************************
*程序用途：用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 CaculateFirstMagnet()                           /*计算初试化后系统磁矩*/
{
    int i,j;
    int iMagnet=0;
    for(i=0;i<L;i++)
    {
        for(j=0;j<L;j++)
   iMagnet += lattice[j];
    }
    magnet = iMagnet;
}

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;
}

C-T曲线：
C-X曲线：

本帖转自52mc论坛