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

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 //设置随机 ...

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


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

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

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

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


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

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

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