isingmodel.h

计算IsingModel的磁化强度

// 调用函数IsingModel2D()即可。
// ---------------------------

// 外部引用
// 随机数发生器，在主程序中初始化
extern CRandom CRandomGenerator;

// 结果记录，在主程序中初始化
extern CLog CResultLog;
// --------------------

// 宏定义
#define SAFE_DELETE_ARRAY(p) { if(p) { delete[] (p);   (p)=NULL; } }
// -----------------------------------------------------------------

// 变量
// 边界条件类型
#define BC_TYPE_FBC 0
#define BC_TYPE_PBC 1

// g_uiLatticeLength 点阵的线度，点阵大小为g_uiLatticeLength的平方
unsigned int g_uiLatticeLength = 0;
unsigned int g_uiLatticeLengthE2 = 0;

// 点阵数组
int** g_ppiLattice = NULL;

// 正方向邻近位置数组。iAdjacentPositive[i]的值为i的正方向邻近点的坐标。
int* g_piAdjacentPositive = NULL;

// 负方向邻近位置数组。
int* g_piAdjacentNegative = NULL;	

// 跃迁几率值。如果大于这个几率就进行反转。
// 实际上，(i,j)与邻近四点相乘求和只有-4,-2,0,2,4几种可能。而只有当2,4的时候才产生随机数，小于0必然跃迁
//	故只有TransitionProbability[2],TransitionProbability[4]才是有效的。
double g_dTransitionProbability[5];

// 计算分析的积累量.
double g_dSumAbsS;				// 序参量S的积累量,注意要积累绝对值
double g_dSumAbsSe2;			// 序参量S的平方积累量,e表示exponential
double g_dSumAbsSe4;			// 序参量S的4次方累积量，为求4阶矩
// --------------------------------

// 函数
// 初始化点阵。这里把所有的点都设成自旋向下
inline void InitLattice( void )
{
	for( unsigned int i = 0; i < g_uiLatticeLength; i++ )
		for( unsigned int j = 0; j < g_uiLatticeLength; j++ )
			g_ppiLattice[i][j] = -1;
}

// 计算跃迁几率。jkt=J/KT。Metropolis函数
inline double CalculateW( double jkt,int sum )		
{
	return exp( (-2.0)*jkt*(double)(sum) );
}

// 蒙特卡罗步。
inline void DoOneMonteCarloStep()
{
	unsigned int i,j;
	int iPrevValue;
	int iSiSjSum;

	// MCS开始
	for( i = 0; i < g_uiLatticeLength; i++ )
		for( j = 0; j < g_uiLatticeLength; j++ )
		{
			iPrevValue = g_ppiLattice[i][j];
			iSiSjSum =(g_ppiLattice[g_piAdjacentPositive[i]][j]
			+ g_ppiLattice[i][g_piAdjacentPositive[j]]
			+ g_ppiLattice[g_piAdjacentNegative[i]][j]
			+ g_ppiLattice[i][g_piAdjacentNegative[j]])*iPrevValue;
			if( iSiSjSum > 0 ) // 根据随机数来决定是否反转
			{
				if( CRandomGenerator.GetRandom() < g_dTransitionProbability[ iSiSjSum ] )
					g_ppiLattice[i][j] = -iPrevValue;
			}
			else	// 能量降低就直接反转
				g_ppiLattice[i][j] = -iPrevValue;
		}
	// MCS结束
}

inline void InitAccVars()
{
	g_dSumAbsS = 0.0;
	g_dSumAbsSe2 = 0.0;
	g_dSumAbsSe4 = 0.0;
}

// 进行采样分析
inline void SampleAnalyse( unsigned long ulCurrStep )
{
	// 计算磁化强度
	unsigned int i,j;
	int iTotalM = 0;
	for( i = 0; i < g_uiLatticeLength; i++ )
		for( j = 0; j < g_uiLatticeLength; j++ )
			iTotalM += g_ppiLattice[i][j];

	// 输出，因为有可能积累的时候改变了符号，所以先输出
	//double dAvgM = ((double)iTotalM)/((double)g_uiLatticeLength)/((double)g_uiLatticeLength);
	//CResultLog.AppendLog( ulCurrStep, dAvgM );

	// 积累量
	if( iTotalM < 0 )
		iTotalM = -iTotalM;	

	double dTemp = (double)iTotalM/(double)g_uiLatticeLengthE2;
	g_dSumAbsS += dTemp;		// 积累一次方
	dTemp *= dTemp;
	g_dSumAbsSe2 += dTemp;		// 积累二次方
	dTemp *= dTemp;
	g_dSumAbsSe4 += dTemp;		// 积累四次方
}

// 完成总步数进行总的计算
inline void GlobalAnalyse( unsigned long ulTotalSamples, double jkt )
{
	// 计算磁化率JX,得出的值其实是H表达式里面的J和磁化率X的乘积
	//CResultLog.WriteDouble( "magnetic susceptibility(multiplied by J)", dJX );
	double dJX = ( (double)(g_dSumAbsSe2) - (double)(g_dSumAbsS*g_dSumAbsS/(double)ulTotalSamples) )/(double)ulTotalSamples *jkt*g_uiLatticeLengthE2;
	CResultLog.WriteDoubleEL( dJX );

	//double d4thM = ((double)(g_ulSumAbsMe4))/((double)(ulTotalSamples));
	//CResultLog.WriteDouble( "fourth moment of intensity of magnetization", d4thM );

	//double d2thM = ((double)(g_ulSumAbsMe2))/((double)(ulTotalSamples));
	//CResultLog.WriteDouble( "second moment of intensity of magnetization", d2thM );
	//double dUL = 1.0-(d4thM/(3.0*d2thM*d2thM));
	//CResultLog.WriteDouble( "accumulated value", dUL );
}

// 收尾工作
void Destroy()
{
	if( g_ppiLattice != NULL )
	{
		for( unsigned int i = 0; i <= g_uiLatticeLength; i++ )
			SAFE_DELETE_ARRAY( g_ppiLattice[i] );
		SAFE_DELETE_ARRAY( g_ppiLattice );
	}

	SAFE_DELETE_ARRAY( g_piAdjacentPositive );
	SAFE_DELETE_ARRAY( g_piAdjacentNegative );
}

// Ising模型运算。

// ulStep 总步数
// ulStepDrop 前面应该舍弃的步数（暂态）
// ulSampleInterval 采样分析（计算磁化强度等）的间隔步数，为了避开相邻两Markov状态之间的关联
void IsingModel2D( unsigned int uiLatticeLength, int iBorderConditionType, unsigned long ulStep, unsigned long ulStepDrop, unsigned long ulSampleInterval, double jkt )
{
	// 假设参数都合法
	g_uiLatticeLength = uiLatticeLength;
	g_uiLatticeLengthE2 = g_uiLatticeLength*g_uiLatticeLength;

	unsigned int i;

	// 点阵数组。
	// ppiLattice[0][0]到ppiLattice[uiLatticeLength-1][uiLatticeLength-1]这个二维点阵是实际有效部分
	// ppiLattice[x][uiLatticeLength]和ppiLattice[uiLatticeLength][x]是用于自由边界条件的,都为0
	g_ppiLattice = new int*[g_uiLatticeLength+1];
	if( g_ppiLattice == NULL )
	{
		CResultLog.WriteString( "Memory Allocation Failed!-g_ppiLattice" );
		Destroy();
		return;
	}
	for( i = 0; i <= g_uiLatticeLength; i++ )
	{
		g_ppiLattice[i] = new int[g_uiLatticeLength+1];
		if( g_ppiLattice[i] == NULL )
		{
			CResultLog.WriteInt( "Memory Allocation Failed!-g_ppiLattice", i );
			Destroy();
			return;
		}
	}
	for( i = 0; i <= g_uiLatticeLength; i++ )
		g_ppiLattice[g_uiLatticeLength][i] = 0;
	for( i = 0; i <= g_uiLatticeLength; i++ )
		g_ppiLattice[i][g_uiLatticeLength] = 0;

	// 初始化邻近位置数据
	g_piAdjacentPositive = new int[g_uiLatticeLength];
	if( g_piAdjacentPositive == NULL )
	{
		CResultLog.WriteString( "Memory Allocation Failed!-g_piAdjacentPositive" );
		Destroy();
		return;
	}
	g_piAdjacentNegative = new int[g_uiLatticeLength];
	if( g_piAdjacentNegative == NULL )
	{
		CResultLog.WriteString( "Memory Allocation Failed!-g_piAdjacentNegative" );
		Destroy();
		return;
	}

	for( i = 0; i < g_uiLatticeLength; i++ )
	{
		g_piAdjacentPositive[i] = i+1;
		g_piAdjacentNegative[i] = i-1;
	}

	// 根据边界条件类型来填边界处的值
	switch( iBorderConditionType )
	{
	case BC_TYPE_FBC:	// g_ppiLattice[x][g_uiLatticeLength] = g_ppiLattice[g_uiLatticeLength][x] = 0
		g_piAdjacentPositive[g_uiLatticeLength-1] = g_uiLatticeLength;
		g_piAdjacentNegative[0] = g_uiLatticeLength;
		break;

	case BC_TYPE_PBC:
		g_piAdjacentPositive[g_uiLatticeLength-1] = 0;
		g_piAdjacentNegative[0] = g_uiLatticeLength-1;
		break;

	default:
		CResultLog.WriteString( "Invalid iBorderConditionType." );
		return;
	}

	// 初始化计算跃迁几率
	g_dTransitionProbability[0] = 0.0;
	g_dTransitionProbability[1] = 0.0;
	g_dTransitionProbability[2] = CalculateW( jkt, 2 );
	g_dTransitionProbability[3] = 0.0;
	g_dTransitionProbability[4] = CalculateW( jkt, 4 );

	// 初始化点阵
	InitLattice();

	// 初始化积累量
	InitAccVars();

	// 计数变量
	unsigned long ulTotalInterval = (unsigned long)( (ulStep-ulStepDrop)/ulSampleInterval );
	unsigned long ulCounter = 0;
	unsigned long ulCurrStep = ulStepDrop;	// 总的当前步数
	unsigned long ulIntervalCounter = 0;

	// 抛弃暂态情况
	for( ulCounter = 0; ulCounter < ulStepDrop; ulCounter++ )
		DoOneMonteCarloStep();

	ulCounter = 0;

	// 主循环开始
	while( 1 )
	{
		DoOneMonteCarloStep();

		ulCounter++;
		if( ulCounter == ulSampleInterval )
		{
			ulCounter = 0;
			ulIntervalCounter++;
			ulCurrStep += ulSampleInterval;
			SampleAnalyse( ulCurrStep );
			if( ulIntervalCounter == ulTotalInterval )
				break;
		}
	}
	// 主循环结束

	// 总结运算
	GlobalAnalyse( ulTotalInterval, jkt );

	// 收尾工作
	Destroy();
}