24小时热门版块排行榜

返回列表

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

[资源] 【转帖】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;

}

下面是用Mathematica处理数据的结果：

[ Last edited by zyj8119 on 2010-9-10 at 17:51 ]

回复此楼

1楼2010-09-10 17:49:48

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

yahoohoo

铁杆木虫 (著名写手)

★ ★ ★ ★ ★ ★ ★
zh1987hs(金币+5):谢谢yahoohoo的精彩回复 2010-09-10 22:54:51
lei0736(金币+2):谢谢 2010-09-11 08:44:08

引用回帖:

Originally posted by zyj8119 at 2010-09-10 17:49:48:
[code]// Copyright(c)2010 Ma Chengzhang <Machengzhang@126.com>
#include
#include
#include //tr1随机数调用
#include //设置随机 ...

这是你自己写的程序吗？很不错啊，:-)。

这是个巨正则系综下的Monte Carlo模拟，因为总的自旋 $ S = \sum_i {s_i} $ 是变化的。

随机数的种子是无所谓精确与否的，代码中的做法只是为了使得不同的模拟使用不同的随机数序列。

建议将程序模块话，比如 Ising MC的一般模块为：

void INIT(); // Set initial configuration

double totalEnergy(); // Calc total energy

double SpinEnergy(int i, int j); // Energy of a given spin

void flip(); // flip a randomly chosen spin

主程序部分的周期性边界条件处理可以更快：

int iRight = i + 1, iLeft = i - 1, jUp = j + 1, jDown = j - 1;
if (iRight = BOXLX) iRight = 0;
if (iLeft = -1) iLeft = BOXLXM1; // BOXLXM1 = BOXLX - 1;
if (jUp = BOXLY) jUp = 0;
if (jDown = -1) jDown = BOXLYM1; // BOXLYM1 = BOXLY - 1;

这样做的好处是将if else 分支结构转变成单独的if，既形式简单又速度更快

能量判断可以简单写成如下代码：

double dE = Enew - Eold;
if (dE <= 0 || rng.gen_open0_open1() < exp(-BETA * dE)) { // 这样的好处是如果dE <= 0 成立那么后面的部分就不用计算
// ACCEPTED
}
else {
// REJECTED
}

关于能量与温度的关系，建议能量用 $ E/NJ $，这样的话在低温区，所有的自旋方向一致，那么构象能为 $ E_{min} = -2 N |J| $。所以最低能量应该为 $ E / NJ = -2 $。而在高温区，自旋向上与向下的概率相等，因为每个自旋近邻的四个自旋和为0，因为总能量也为0。根据这两点，模拟得到的结果是有问题的。

[ Last edited by yahoohoo on 2010-9-10 at 21:13 ]

赞一下(2人)

回复此楼

2楼2010-09-10 21:05:59

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

zyj8119

木虫 (著名写手)

引用回帖:

Originally posted by yahoohoo at 2010-09-10 21:05:59:

这是你自己写的程序吗？很不错啊，:-)。

这是个巨正则系综下的Monte Carlo模拟，因为总的自旋 $ S = \sum_i {s_i} $ 是变化的。

随机数的种子是无所谓精确与否的，代码中的做法只是为了使得不同的模拟使 ...

谢谢建议！！

回复此楼

3楼2010-09-10 22:56:11

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖