ising模型c语言程序,【转帖】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<

//判断是否反转自旋

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<

//cout<

///

cout<

char c;

cin>>c;

}