练习 3. Ising 模型和 Metropolis 算法¶
这 不是 一个传统的 OJ 题目，而更倾向于一个计算物理实践，一步一步地完成编程步骤而达到最后的目标。如果你希望了解一个从物理模型上构建的 C 程序项目可以怎样编写，可以跟着实践。
步骤 1. 关于 Ising 模型¶
伊辛模型（Ising model），是一个以物理学家恩斯特·伊辛为名的数学模型，用于描述物质的铁磁性。该模型中包含了可以用来描述单个原子磁矩的参数 \(\sigma _{i}\) ，其值只能为\(+1\)或\(-1\)，分别代表自旋向上或向下，这些磁矩通常会按照某种规则排列，形成晶格，并且在模型中会引入特定交互作用的参数，使得相邻的自旋互相影响。虽然该模型相对于物理现实是一个相当简化的模型，但它却和铁磁性物质一样会产生相变。事实上，一个二维的方晶格伊辛模型是已知最简单而会产生相变的物理系统。
对于两个相邻的晶格点 \(i,j\)，我们可以引入一个交互作用参数 \(J_{ij}\)，此外，我们可以假设每个自旋 \(j\) 都和外加的磁场 \(h_j\) 作用。则整个系统的哈密顿量（可以理解为总能量）可写成:
其中 \(<ij>\) 代表晶格点 i 和晶格点 j 是相邻的晶格点。因此哈密顿量的第一项为对每一对相邻晶格点的总和（每一对只算一次），代表所有自旋之间交互作用的能量，而第二项则是磁场和自旋交互作用的能量。µ 是晶格点磁矩的值，值得注意的是，电子的磁矩和他的自旋方向相反，所以哈密顿量的第二项应该要是正号比较合理，但在习惯上，还是会令第二项为负号。
一个常见的简化是假设没有外加的磁场作用在伊辛模型上，也就是说，对于所有的 \(j∈Λ，h_j = 0\)。利用这项简化，其哈密顿量可以写成:
另一个常见的简化是假设所有相邻晶格点的交互作用都是相等的，因此对于所有相邻的 \(i\) 和 \(j\)，可以假设 \(J_{ij} = J\)，而其哈密顿量可以写为:
步骤 2. 基本数据结构¶
我们希望进行二维 Ising 模型的计算模拟，因此需要建立一个二维数组来记录 \(\sigma_{(x,y)}\) 的值。由于每个晶格点只能取向上（\(+1\)）或者向下（\(-1\)）两种状态，因此一个 int8_t 已经足够了。
#define WIDTH 100
#define HEIGHT 100
uint8_t sigma[WIDTH][HEIGHT];
步骤 3. 体系 Hamilton 量的计算¶
在下面的演算过程我们会采用周期性边界条件使的每个晶格点的相邻数 \(c\) 都相等。这里需要注意的是，交互项 \(J\sigma_i\sigma_j\) 对于每对 \((i, j)\) 只需要计算一次。
#define J 1.0
double hamilton()
{
double h = 0; // Hamilton 求和
for (int i = 0; i < WIDTH; ++i)
{
for (int j = 0; j < HEIGHT; ++j)
{
int nexti = (i + 1) % 100; // 周期边界条件
int nextj = (j = 1) % 100;
h += -J * (sigma[i][j] * sigma[i + 1][j]);
h += -J * (sigma[i][j] * sigma[i][j + 1]);
}
}
return h;
}
步骤 4. Metropolis 算法¶
Metropolis 算法是在数值模拟伊辛模型时最常用的一种蒙特卡洛方法。这种方法首先要选定一个选择几率 \(g(\mu, \nu)\)，代表系统在状态 \(\mu\) 下，在所有可能状态中选到状态 \(\nu\) 的几率。另外还需定义出一个接受几率 \(A(\mu, \nu)\) ，也就是说当系统在状态 \(\mu\) 下，接受系统跳到态 \(\nu\) 的几率。如此的设计是为了让系统达到细致平衡。如果状态 \(\nu\) 被接受了，则整个系统就会跳到状态 \(\nu\)，并且选择和决定下一个要跳到的状态。如果状态 \(\nu\) 被没被接受，系统则留在状态 \(\mu\)，一样重新选择下一个要跳到的状态。这样的步骤一直重复直到某些条件达成为止，譬如说整个伊辛模型完全被磁化，也就是所有的自旋都只到同一个方向。另外在实行这种算法，有一点需要注意的是需要选到适当 \(g(\mu, \nu)\) 以保证整个过程的遍历性(ergodicity)。
在热平衡时，整个系统的能量只会有小幅度的扰动，这点促成了在演算时采用单一自旋反转法进行计算，也就是说每次系统转换其状态时，只改变其中一个自旋的方向。对自旋数很多的一伊辛模型来说，系统在不同的状态之间跳跃时，其能量改变的幅度都很小。事实上，对于每个晶格点都和 \(c\) 个晶格点相邻的模型来说，每次能量改变的最大幅度为 \(2cJ\)。此外，采用这种单一自旋反转法可以保证演算过程的遍历性，因为任意一个状态都可以借由逐次的反转相异的自旋，而变成任意其他状态。
因为共有 \(L\) 个晶格点，而且单一自旋反转法是唯一系统可以从一个状态跳到另一个状态的方法，所对于任何一个状态 \(\mu\)，有 \(L\) 个新状态 \(\nu\) 它可以跳去。可以假设 \(\mu\) 跳到这 \(L\) 个新状态的几率是相等的，因此 \(g(\mu, \nu) = 1/L\)。 为了满足细致平衡，所以下面的等式要成立:
因此我们希望接受几率满足:
如果 \(H_\nu > H_\mu\) 则 \(A(\nu, \mu) > A(\mu, \nu)\)，因此 Metropolis 令 \(A(\mu, \nu)\) 和 \(A(\nu, \mu)\) 中较大的为 \(1\)。按照这样的定可以得到 \(A(\mu, \nu)\) 的值为:
最后，整个算法最基本的形式为:
- 用 \(g(\mu, \nu)\) 选出一个自旋，并且计算所有和其自旋相关的能量贡献。(也就是所有和它相邻自旋交互作用的能量贡献。)
- 反转其自旋，再计算一次所有和它相关的能量贡献。
- 如果其能量贡献下降，保持自旋反转。
- 如果贡献能量上升，则令自旋有 \(e^{-\beta (H_\nu-H_\mu)}\) 的几率保持反转。
- 重复步骤一。
整个系统能量的改变量 \(H_\nu−H_\mu\) 只和反转的自旋和和它相邻的自旋有关，所以只要相邻的自旋数不要太多，计算的速度算是相当快的。而整个系统会逐渐的趋近于某个平衡。
步骤 5. 计算一个自旋的 Hamilton 量贡献¶
上面说到，能量的改编量其实只与相邻的自旋状态有关，而我们只关心 \(H_\mu\) 和 \(H_\nu\) 的差距，不关心它们的绝对大小。因此，可以只计算一个自旋的贡献。
double hamilton_contribution(int i, int j, int8_t state)
{
// 应用周期条件
int iminus = i == 0 ? WIDTH - 1 : i - 1;
int iplus = i == WIDTH - 1 ? 0 : i + 1;
int jminus = j == 0 ? HEIGHT - 1 : j - 1;
int jplus = j == HEIGHT - 1 ? 0 : j + 1;
double hcontrib = -J * state * (
sigma[iminus][j] +
sigma[iplus][j] +
sigma[i][jminus] +
sigma[i][jplus]
);
return hcontrib;
}
步骤 6. 实现随机数的选取¶
从上面的分析可以知道我们需要生成 \([0, \texttt{WIDTH})\) 和 \([0, \texttt{HEIGHT})\) 的随机整数，以及 \([0, 1)\) 的随机小数。因此，需要将系统的 rand() 函数的值进行一定的修整。注意，直接使用类似 rand() % WIDTH 这样的表达式是不正确的，因为这导致随机数的不均匀。
double rand_01()
{
return rand() / (double)(RAND_MAX + 1);
}
int rand_range(int low_inclusive, int high_exclusive)
{
return (int)(rand_01() * (high_exclusive - low_inclusive)) + low_inclusive;
}
步骤 7. 实现 Metropolis 算法的一步¶
按照基本思路，Metropolis 算法是迭代性质的，因此，先来实现它的一步。
void metropolis_step(double beta)
{
// 选取随机的一个咨自旋
int randi = rand_range(0, WIDTH);
int randj = rand_range(0, HEIGHT);
// 计算这个自旋现在的状态带来的 Hamilton 量
double h_current = hamilton_contribution(randi, randj, sigma[i][j]);
// 计算这个自旋翻转带来的 Hamilton 量
double h_inverse = hamilton_contribution(randi, randj, -sigma[i][j]);
if (h_inverse < h_current)
{
// 如果翻转能量更小，那么保持翻转
sigma[i][j] = -sigma[i][j];
}
else
{
// 以一定的概率接受能量提高的翻转
if(rand_01() < pow(M_E, -beta * (h_inverse - h_current)))
{
sigma[i][j] = -sigma[i][j];
}
}
}
步骤 8. 完成整个程序并加入可视化¶
#include <stdio.h>
#include <math.h>
#include <stdint.h>
#include <stdlib.h>
#include <time.h>
#define WIDTH 100
#define HEIGHT 100
uint8_t sigma[WIDTH][HEIGHT];
#define J 1.0
double rand_01()
{
return rand() / (double)(RAND_MAX + 1);
}
int rand_range(int low_inclusive, int high_exclusive)
{
return (int)(rand_01() * (high_exclusive - low_inclusive)) + low_inclusive;
}
double hamilton_contribution(int i, int j, int8_t state)
{
int iminus = i == 0 ? WIDTH - 1 : i - 1;
int iplus = i == WIDTH - 1 ? 0 : i + 1;
int jminus = j == 0 ? HEIGHT - 1 : j - 1;
int jplus = j == HEIGHT - 1 ? 0 : j + 1;
return -J * state * (
sigma[iminus][j] +
sigma[iplus][j] +
sigma[i][jminus] +
sigma[i][jplus]
);
}
void metropolis_step(double beta)
{
int i = rand_range(0, WIDTH);
int j = rand_range(0, HEIGHT);
double h_current = hamilton_contribution(i, j, sigma[i][j]);
double h_inverse = hamilton_contribution(i, j, -sigma[i][j]);
if (h_inverse < h_current || rand_01() < exp(-beta * (h_inverse - h_current)))
{
sigma[i][j] = -sigma[i][j];
}
}
void visualize()
{
for (int i = 0; i < WIDTH; ++i)
{
for (int j = 0; j < HEIGHT; ++j)
{
printf("%s", sigma[i][j] == 1 ? "##" : " ");
}
printf("\n");
}
printf("---------\n");
}
int main(void)
{
srand((unsigned int)time(NULL)); // Seed RNG
double beta = 0.1; // Change as needed
// Initialize spins to ±1
for (int i = 0; i < WIDTH; ++i)
{
for (int j = 0; j < HEIGHT; ++j)
{
sigma[i][j] = rand_01() > 0.5 ? 1 : -1;
}
}
// Run Monte Carlo steps
for (int step = 0; step < 10000; ++step)
{
metropolis_step(beta);
if (step % 100 == 0) // visualize every 100 steps
{
visualize();
}
}
return 0;
}
结束语. 观察计算结果¶
- 由于 \(\beta = 1/(kT)\)，调节 \(\beta\) 的值可以调节“温度”，直接决定接受能量提高的翻转的概率。试一试，改变 \(\beta\) 的值，观察一定步数之后体系的状态。提示：\(\beta\) 的临界值约为 \(0.44\)。
- 试一试给体系增加外部交互项，也就是在 \(H\) 中加入一项，得到
在这种情况下，改变 \(h\) 的值，观察结果的变化。
- 现在，初始情况的 \(\sigma_{ij}\) 是随机的。如果初始化为某些具体的内容呢？例如，只有中心点为向上自旋，其余为向下，会发生什么？