1. 程式人生 > 其它 >C語言實現二維伊辛模型的蒙特卡羅方法模擬

C語言實現二維伊辛模型的蒙特卡羅方法模擬

技術標籤:實驗筆記c語言

伊辛模型

簡介

伊辛模型(Ising model)是一類描述物質相變的隨機過程(stochastic process)模型。物質經過相變,要出現新的結構和物性。發生相變的系統一般是在分子之間有較強相互作用的系統,又稱合作系統。

起源

伊辛模型由德國物理學家威廉·楞次(Wilhelm Lenz)在1920年提出以描述鐵磁性物質的內部的原子自旋狀態及其與巨集觀磁矩的關係。1924年,楞次的學生Ernst Ising求解了不包含相變的一維伊辛模型。

闡述

伊辛模型假定分子排布為週期性點陣。在一維、二維、三維空間分別等間距排布、正方形排布、正方體排布;每個分子都有自旋(固有磁矩);系統中可能存在外磁場;系統的Hamilton量(可以看做能量)由外場能和相互作用能共同構成:

H = − J ∑ < i , j > σ i σ j − μ H ∑ k σ k H = -J \sum_{<i,j>}{\sigma_i \sigma_j}-\mu H \sum_{k}{\sigma_k} H=J<i,j>σiσjμHkσk
第一項表示相互作用:相互作用的是兩個最鄰近的分子,並對所有這樣的分子對進行求和。第二項表示外場:μ為玻爾磁子。
Ising模型的特色在於考慮到了分子間的相互作用;但它只考慮最鄰近的分子間的相互作用(比如在二維情形中,則不考慮位於對角線的分子間的相互作用)。

3x3Ising模型格點相互作用
如上圖,為3x3規模的Ising模型,圖中的連線表示系統格點的所有相互作用,符合週期性邊界條件。

凝練《熱力學與統計物理》中關於伊辛模型的介紹,設沒有外磁場,伊辛模型系統的哈密頓量為:
H = − J ∑ < i , j > σ i σ j H = -J \sum_{<i,j>}{\sigma_i \sigma_j} H=J<i,j>σiσj
其中J為相互作用常數,J>0表示系統自旋之間為鐵磁相互作用,σ為自旋方向,取值為-1或1。

蒙特卡羅方法

簡介

蒙特卡羅方法又稱統計模擬法、隨機抽樣技術,是一種隨機模擬方法,以概率和統計理論方法為基礎的一種計算方法,是使用隨機數(或更常見的偽隨機數)來解決很多計算問題的方法。將所求解的問題同一定的概率模型相聯絡,用電子計算機實現統計模擬或抽樣,以獲得問題的近似解。為象徵性地表明這一方法的概率統計特徵,故借用賭城蒙特卡羅命名。

起源

蒙特卡羅方法於20世紀40年代美國在第二次世界大戰中研製原子彈的“曼哈頓計劃”計劃的成員S.M.烏拉姆和J.馮·諾伊曼首先提出。數學家馮·諾伊曼用馳名世界的賭城—摩納哥的Monte Carlo—來命名這種方法,為它蒙上了一層神祕色彩。在這之前,蒙特卡羅方法就已經存在。1777年,法國Buffon提出用投針實驗的方法求圓周率。這被認為是蒙特卡羅方法的起源。

闡釋

蒙特卡羅方法的基本思想是當所求解問題是某種隨機事件出現的概率,或者是某個隨機變數的期望值時,通過某種“實驗”的方法,以這種事件出現的頻率估計這一隨機事件的概率,或者得到這個隨機變數的某些數字特徵,並將其作為問題的解。
蒙特卡羅方法可以用來模擬系統的物理性質,如自發磁化強度,其定義為:
M = ∑ i σ i N M=\frac{\sum_i{\sigma_i}}{N} M=Niσi
本實驗將利用蒙特卡羅方法來模擬Ising模型表示的物理系統的隨機變化過程。模擬過程的流程圖如下:

選定自旋 翻轉該自旋 ΔE>0? exp(-ΔE/kT)<R? 接受改變

選定系統任意初態,按照上述方法實現自旋翻轉,並遍歷所有自旋,這一過程稱為一個蒙特卡羅步。其中,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);
}