24小时热门版块排行榜

返回列表

【奖励】本帖被评价4次，作者zyj8119增加金币 3 个

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

zyj8119

木虫 (著名写手)

[资源] 【转帖】2维无外场Ising模型的Monte Carlo模拟的C++实现

CODE:

// Copyright(c)2010 Ma Chengzhang 

#include 

#include 

#include  //tr1随机数调用

#include  //设置随机数种子用

#include 

using namespace std; //cout,endl等简单语句必须

using namespace std::tr1; //tr1随机数调用

inline unsigned __int64 ReadTSC()//为取到精确的随机数种子，嵌入汇编程序取系统时钟

{

 __asm _emit 0x0F

  __asm _emit 0x31

}

double totalE(vector >& SpinGrid,int GridL,double JBound)

{

 int nn_count = 0;//最近邻计数

 // 每个格点只考虑右键和上键

 for(int i=0; i
 {

  for(int j=0; j
  {

   if(SpinGrid[i][j]==SpinGrid[i+1][j]) // 考虑右键

    nn_count++;

   else

    nn_count--;

   if(SpinGrid[i][j]==SpinGrid[i][j+1]) // 考虑上键

    nn_count++;

   else

    nn_count--;

  }

 }

 //周期性边界情况

 for(int i=0; i
 {

  if(SpinGrid[i][GridL-1]==SpinGrid[i+1][GridL-1]) // 右键

   nn_count++;

  else

   nn_count--;

  if(SpinGrid[i][GridL-1]==SpinGrid[i][0]) // 上键

   nn_count++;

  else

   nn_count--;

 }

 for(int j=0; j
 {

  if(SpinGrid[GridL-1][j]==SpinGrid[0][j]) // 右键

   nn_count++;

  else

   nn_count--;

  if(SpinGrid[GridL-1][j]==SpinGrid[GridL-1][j+1]) // 上键

   nn_count++;

  else

   nn_count--;

 }

 //右上顶点

 if(SpinGrid[GridL-1][GridL-1]==SpinGrid[0][GridL-1]) // 右键

  nn_count++;

 else

  nn_count--;

 if(SpinGrid[GridL-1][GridL-1]==SpinGrid[GridL-1][0]) // 上键

  nn_count++;

 else

  nn_count--;

 return -JBound*nn_count;

}

void main()

{

 //物理常数定义

 double kB=1.0; //Boltzmann常数

 double JBound=0.5; //交换作用系数

 int GridL=20; //网格大小，长宽一样，即 GridL x GridL 网格

 mt19937 RandEngine; //Mersenne Twister算法随机数生成器，核心引擎类

 unsigned long RandSeed = (unsigned long)ReadTSC(); //用 CPU 时钟作为随机数种子

 RandEngine.seed(RandSeed);

 uniform_int RandInt(0, GridL-1);//产生0到GridL-1之间的随机整数,

 // 含0和GridL-1，即[0,GridL-1]而不是(0,GridL-1)

 uniform_real RandFloat(0,1); //产生0到1之间的随机数

 double MaxTemp=5.00; //要模拟的最高温度，温度以J/kB为单位，kB是波尔兹曼常数

 double MinTemp=0.1; //要模拟的最低温度，温度以J/kB为单位，kB是波尔兹曼常数

 unsigned int BeforeEquilSteps=10000; //在达到平衡前要演化的步数

 unsigned int TotalSteps=BeforeEquilSteps+10000; //每温度要演化的总步数，包含前面的BeforeEquilSteps

 int DiscardedSteps=50; //达到平衡态后，为实现无偏统计抽样，每DiscardedSteps演化步数才统计一次

 double TempInterval=0.1; //扫描温度间隔，温度以J/kB为单位，kB是波尔兹曼常数

 double beta; //当前beta=1/(kB*T)

 vector< vector > SpinGrid(GridL, vector(GridL,-1)); //网格定义，格点初始值都取-1

 //统计量定义

 double TotalEnergy=0; //网格总能量

 double persiteEnergy=0; //每格点能量

 double persiteAveEnergy=0; //每格点平均能量

 ofstream fileout("Ising2D.dat");

 //Monte Carlo 部分

 for(double CurrentTemp=MinTemp;CurrentTemp
 {

  int count=0; //此变量记录达到平衡后演化步数，每当count达到DiscardedSteps时归零

  int Outcount=0; //记录累计统计次数

  int left,right,bottom,top; //辅助变量

  beta=1/(kB*CurrentTemp);

  for(unsigned int mcStep=0;mcStep
  {

   //随机选择格点反转自旋，随机数算法取自：C++ tr1 Mersenne Twister算法

   int i=RandInt(RandEngine); //第i列

   int j=RandInt(RandEngine); //第j行

   //采用周期性边界条件

   if(i!=0)

    left = SpinGrid[i-1][j];

   else

    left = SpinGrid[GridL-1][j];

   if(i!=GridL-1)

    right = SpinGrid[i+1][j];

   else

    right = SpinGrid[0][j];

   if(j!=0)

    bottom = SpinGrid[i][j-1];

   else

    bottom = SpinGrid[i][GridL-1];

   if(j!=GridL-1)

    top = SpinGrid[i][j+1];

   else

    top = SpinGrid[i][0];

   int spinSum=top+bottom+left+right;

   //double oldEnergy=-JBound*SpinGrid[i][j]*spinSum;

   //double newEnergy=-oldEnergy;

   double EDiff=2*JBound*SpinGrid[i][j]*spinSum;//EDiff=newEnergy-oldEnergy

   //cout<<"EDiff="<

   //判断是否反转自旋

   if( EDiff<0)

    SpinGrid[i][j]=-SpinGrid[i][j];

   else

   {

    if(RandFloat(RandEngine)
     SpinGrid[i][j]=-SpinGrid[i][j];

   }

   //如果mcStep达到了指定步数，就认为达到了平衡态，可以进行状态分析了

   //但是只是每隔DiscardedSteps步数才分析一次

   if(mcStep>BeforeEquilSteps)

   {

    count=count+1;

    if(count==DiscardedSteps)

    {

     count=0; //count归零，重新开始计数

     TotalEnergy=totalE(SpinGrid,GridL,JBound);

     persiteEnergy=TotalEnergy/(GridL*GridL);  

     persiteAveEnergy=(persiteAveEnergy*Outcount+persiteEnergy)/(Outcount+1);

     Outcount++;

    }

   }

  }

  //输出数据到数据文件

  fileout<< CurrentTemp <<" "<

  //统计量归零

  persiteAveEnergy=0;

  Outcount=0;

 }

 //Monte Carlo 部分结束

 fileout.close();

 ////输出测试

 ////测试随机整数

 //for(int i = 0; i < 50; ++i) //产生均匀分布的在0到19之间的50个整数随机数

 //{

 // std::cout << RandInt(RandEngine) << endl;

 //}

 //cout<<"*************"<
 //for(int i = 0; i < 50; ++i) //产生均匀分布的在0到1之间的50个随机数

 //{

 // std::cout << RandFloat(RandEngine) << endl;

 //}

 ////测试随机整数毕

 //cout<
 //cout<<"TotalSteps = "<
 //cout<<"RandSeed = "<
 ///////////////////////////////////////////////////////////////////////

 cout<
 char c;

 cin>>c;

}