isingb.cpp

模拟计算磁滞回线

#include "fstream.h"
#include "math.h"
#include "stdlib.h"
#include "iostream.h"

double JKT =0.30;
//完成测试
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 ubkt)
{
	//de能量差
	int i,m=0;
	double x,y,p,h,t,de;
	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=(JKT*m+ubkt)*(-2)*a[xx*n+yy];//+ub?-ub?  de mean -de/kt
		p=exp(de);//
		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 ubkt,int total)
{
	for(int i=0;i<total;i++)
	{
		MCS(a,n,ubkt);
	}
}
//抽样count次 每隔10 MCS 抽样一次 返回值为count次抽样的平均值
double sample(int a[],int n,double ubkt,int count)
{	
	double smp=0.0;
	for(int i=0;i<count;i++)
	{
		MCStime(a,n,ubkt,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,ubkt=1.0;

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

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