ising.cpp

Ising模型

#include "fstream.h"
#include "math.h"
#include "stdlib.h"
#include "iostream.h"
//完成测试
void initialize(int a[],int n)
{	
	int i,j;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			int long x=rand();
			if(x<RAND_MAX/2)
				a[i*n+j]=1;
			else 
				a[i*n+j]=-1;
		}
	}
}
/*
//MCS
void MCS(int a[],int n,double jkt)
{
	//de能量差
	int i,j,m=0,de;
	double p,h;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			m =0;
			m+=a[i*n+(j+1)%n]; //右
			m+=a[i*n+(j+n-1)%n]; //左
			m+=a[((i+1)%n)*n+j]; //下
			m+=a[((i+n-1)%n)*n+j]; //上
			de=m*2*a[i*n+j]; 
		//	temp=max(0,de);// for -if
		//	p=exp(-1*temp*jkt);
			p=exp(-1*de*jkt);
			if(p<1)
			{
				h=rand()*1.0/RAND_MAX;
				if(h<p) a[i*n+j]=a[i*n+j]*(-1);
			}
			else a[i*n+j]=a[i*n+j]*(-1);
		}
	}
}
*/
//MCS
void MCS(int a[],int n,double jkt)
{
	//de能量差
	int i,m=0,de;
	double x,y,p,h,t;
	t=1.0/n; // t步长
	for(i=0;i<n*n;i++)
	{	
		
		x=rand()*1.0/RAND_MAX;
		y=rand()*1.0/RAND_MAX;
		int xx=int(x/t);
		int yy=int(y/t);
		m =0;
		m+=a[xx*n+(yy+1)%n]; //右
		m+=a[xx*n+(yy+n-1)%n]; //左
		m+=a[((xx+1)%n)*n+yy]; //下
		m+=a[((xx+n-1)%n)*n+yy]; //上
		de=m*2*a[xx*n+yy];
		p=exp(-1*de*jkt);
		if(p<1)
		{
			h=rand()*1.0/RAND_MAX;
			if(h<p) a[xx*n+yy]=a[xx*n+yy]*(-1);
		}
		else a[xx*n+yy]=a[xx*n+yy]*(-1);
	}

}
	
//total次MCS
void MCStime(int a[],int n,double jkt,int total)
{
	for(int i=0;i<total;i++)
	{
		MCS(a,n,jkt);
	}
}
//抽样count次 每隔10 MCS 抽样一次 返回值为count次抽样的平均值
double sample(int a[],int n,double jkt,int count)
{	
	double smp=0.0;
	for(int i=0;i<count;i++)
	{
		MCStime(a,n,jkt,10);
		for(int j=0;j<n;j++)
		{
			for(int k=0;k<n;k++)
			{
				smp+=a[j*n+k];
			}
		}
	}
	smp=smp/(n*n)/count;
	return smp;
}
void main()
{
	int total=2000,n=40,count=100;
	double s;

	int * p;
	p=new int[n*n];
	//初始化
	initialize(p,n);

	ofstream outfile("Ising.dat");
	if(! outfile)
	{
		cout<<"can't open the file";
		exit(1);
	}
	
	for(double jkt=0.40;jkt<0.50;jkt+=0.002)
	{
		//多次MCS 达到平衡态
		MCStime(p,n,jkt,total);
		//平衡后抽样，s为抽样平均值
		s=sample(p,n,jkt,count);
		//输出当前jkt下的<s>
		outfile<<jkt<<'\t'<<s<<endl;
	}
	outfile.close();
	delete p;
}