物理：C语言实现二维伊辛模型的蒙特卡罗方法模拟

伊辛模型

简介

伊辛模型（Ising model）是一类描述物质相变的随机过程（stochastic process）模型。物质经过相变，要出现新的结构和物性。发生相变的系统一般是在分子之间有较强相互作用的系统，又称合作系统。

起源

伊辛模型由德国物理学家威廉·楞次（Wilhelm Lenz）在1920年提出以描述铁磁性物质的内部的原子自旋状态及其与宏观磁矩的关系。1924年，楞次的学生Ernst Ising求解了不包含相变的一维伊辛模型。

阐述

伊辛模型假定分子排布为周期性点阵。在一维、二维、三维空间分别等间距排布、正方形排布、正方体排布；每个分子都有自旋（固有磁矩）；系统中可能存在外磁场；系统的Hamilton量（可以看做能量）由外场能和相互作用能共同构成：
$$E=-J\sum _{<i,j>}{\sigma }_i{\sigma }_j-\mu H\sum _k{\sigma }_k$$
第一项表示相互作用：其中J为相互作用常数，J>0表示系统自旋之间为铁磁相互作用，σ为自旋方向，取值为-1或1。参与相互作用的是两个最邻近的分子，并对所有这样的分子对进行求和。第二项表示外场：μ为玻尔磁子。
系统的磁化强度定为格点的平均自旋方向：
$$M=\frac{{\sum }_i{\sigma }_i}{N}$$
其中N是系统的格点数。
Ising模型的特色在于考虑到了分子间的相互作用；但它只考虑最邻近的分子间的相互作用（比如在二维情形中，则不考虑位于对角线的分子间的相互作用）。

如上图，为3x3规模的Ising模型，图中的连线表示两个格点间的存在相互作用，要符合周期性边界条件。
凝练《热力学与统计物理》中关于伊辛模型的介绍，设没有外磁场，则伊辛模型系统的哈密顿量为：
$$H=-J\sum _{<i,j>}{\sigma }_i{\sigma }_j$$
没有外磁场影响下的磁化强度为自发磁化强度。

蒙特卡罗方法

简介

蒙特卡罗方法又称统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。

起源

蒙特卡罗方法于20世纪40年代美国在第二次世界大战中研制原子弹的“曼哈顿计划”计划的成员S.M.乌拉姆和J.冯·诺伊曼首先提出。数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种方法，为它蒙上了一层神秘色彩。在这之前，蒙特卡罗方法就已经存在。1777年，法国Buffon提出用投针实验的方法求圆周率。这被认为是蒙特卡罗方法的起源。

阐释

蒙特卡罗方法的基本思想是当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。
蒙特卡罗方法可以用来模拟系统的物理性质。
本实验将利用蒙特卡罗方法来模拟Ising模型表示的物理系统的随机变化过程。模拟过程的流程图如下：

选定系统任意初态，按照上述方法实现自旋翻转，并遍历所有自旋，这一过程称为一个蒙特卡罗步。其中，R为0~1间的一个实数，来决定翻转概率为exp(-ΔE/kT)的自旋是否翻转。重复该过程可使系统达到（能量）比较稳定的状态。

C语言代码

//预处理部分
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#pragma warning(disable:4996)
//忽略VS开发环境下的安全性错误提示C4996

//生成一个正方形矩阵，元素随机赋值-1或1
int** createAMatrix(int theOrder) {
	int** aMatrix = (int**)malloc(sizeof(int*) * theOrder);
	int i, j;
	for (i = 0; i < theOrder; i++) {
		aMatrix[i] = (int*)malloc(sizeof(int) * theOrder);
		for (j = 0; j < theOrder; j++)aMatrix[i][j] = (rand() % 2) * 2 - 1;
	}
	return aMatrix;
}

//计算正方形矩阵所表示的ising模型的平均磁化强度
//Map为矩阵的指针
//Size为矩阵的边长
double theM(int** aMatrix, int theOrder) {
	int i, j, SUM = 0;
	for (i = 0; i < theOrder; i++)for (j = 0; j < theOrder; j++)SUM += aMatrix[i][j];
	return SUM * 1.0 / theOrder / theOrder;
}

//对正方形矩阵所表示的ising模型中的每个点按顺序进行蒙特卡罗模拟优化
//Map为矩阵的指针
//Size为矩阵的边长
//T为ising模型的温度
//J为ising模型中各点的相互作用系数
void monteCarlo(int** aMatrix, int theOrder, double T, double J) {
	int i, j, E, dE;
	double R;
	for (i = 0; i < theOrder; i++)for (j = 0; j < theOrder; j++) {
		dE = 2 * J * aMatrix[i][j] * (aMatrix[(i + theOrder - 1) % theOrder][j] + aMatrix[i][(j + theOrder + 1) % theOrder] + aMatrix[(i + theOrder + 1) % theOrder][j] + aMatrix[i][(j + theOrder - 1) % theOrder]);
		//+ theOrder 及 % theOrder 是为了在边界进行周期性循环
		R = rand() * 1.0 / RAND_MAX;
		if ((dE <= 0) || (exp(-dE * 1.0 / T) >= R))aMatrix[i][j] = -aMatrix[i][j];
	}
}

//清除正方形矩阵（占用的内存）
//Map为矩阵的指针
//Size为矩阵的边长
void clearMatrix(int** aMatrix, int theOrder) {
	int i;
	for (i = 0; i < theOrder; i++)free(aMatrix[i]);
	free(aMatrix);
}

//模拟磁化强度随温度变化的主函数，数据显示到屏幕并保存到M-T_Data.txt文件中
void main() {
	int i, j, theOrder = 20, N = 1000;
	double T, SUM, J = 1.0, L = 0.01;
	int** aMatrix;
	FILE* fp = NULL;
	fp = fopen("Data.txt", "w");
	for (T = 0.1; T <= 5; T += L) {
		SUM = 0;
		aMatrix = createAMatrix(theOrder);
		for (i = N; i > 0; i--) {
			monteCarlo(aMatrix, theOrder, T, J);
			if (i <= 500)SUM += fabs(theM(aMatrix, theOrder));
			//对后500个蒙特卡罗步的系统磁化强度取平均
		}
		clearMatrix(aMatrix, theOrder);
		printf("%.2lf\t%.5lf\n", T, SUM / 500);
		fprintf(fp, "%.2lf\t%.5lf\n", T, SUM / 500);
	}
	fclose(fp);
}