声明：旨在自学一些物理概念和计算方法后整理的笔记和心得，本文基于Kotze先生的Introduction to Monte Carlo methods for an Ising Model of a Ferromagnet 文章做了机械无脑式的翻译，不做任何商业用途，仅供学习交流之用。

原文链接https://arxiv.org/abs/0803.0217v1​arxiv.org

Kotze的这篇文章主要分三个部分：

一、理论和背景介绍

二、计算结果展示

三、程序代码

首先，

1 理论

1.1 引言

文章讨论是对使用蒙特卡罗模拟来评估铁磁体的可观察性的一种有用方法的介绍。主要背景是关于这种随机方法的相关性和有效性，特别是Metropolitan - Hastings算法的适用性。重要的是，本文强调了自发磁化的潜在破坏性影响，并研究了避免这种影响的方法。

引入了一种伊辛模型，用于研究二维铁磁体在不同温度下的磁化强度和能量的性质。计算了观测值，并对临界温度下的相变作了说明和评价。最后通过有限尺度尺度分析确定临界指数，并利用不同格点尺寸的累积量比计算出居里温度。模拟结果与精确计算结果进行了比较，验证了数值计算过程的正确性。所使用的代码的副本，用c++编写的，在通用公共许可证下可免费使用和修改。

1.2 背景

在大多数普通材料中，相关的原子磁偶极子的方向是随机的。实际上，这种非特异性分布不存在整体的宏观磁矩。然而，在某些情况下，如铁，磁矩产生的结果，首选排列的原子自旋。

这一现象基于两个基本原理，即能量最小化和熵最大化。这些是相互竞争的原则，在缓和总体影响方面很重要。温度是这些对立元素之间的中介，并最终决定哪一个更占优势。

能量最小化和熵最大化的相对重要性在本质上是由一个特定的概率决定的

...(1)

如图1所示，称为吉布斯分布。图1，该图展示了变化能量(E)和温度(T)的玻尔兹曼概率分布。

1.3 模型

理解这个理论的一个关键因素是自旋和相关的磁矩。知道了自旋是一种量子力学现象，就很容易预测，对这个问题的完整而彻底的阐述将需要自旋和角动量的量子规则。这些因素被证明是不必要的复杂因素。

因此，我们引入一个模型来获得有用的结果。模型的核心思想是将问题的复杂性简化到数学上易于处理的程度，同时保留系统的基本物理性质。伊辛模型非常有效地做到了这一点，甚至允许一个良好的概念性理解。图2，一个伊辛模型的二维格图。向上和向下的箭头分别表示一个正自旋和一个负自旋。

伊辛模型考虑了二维的问题，并将偶极子自旋放置在规则的格点上，同时限制其自旋轴为上(+y)或下(-y)。晶格构型是正方形，尺寸为L，自旋总数等于

。在最简单的形式中，偶极子之间的相互作用范围被限制在最近邻的位置。这就产生了一个特定自旋位置的哈密顿量，Hi

...(2)

最近邻之间的耦合常数用 J 表示，

和

分别为最近邻自旋。模型中相互作用的性质全部包含在相互作用耦合常数J的符号中。如果J是正的，就意味着该材料具有铁磁性(平行排列)，而负号则意味着该材料是反铁磁性的(倾向于反平行排列)。在我们的讨论中，J 将被取为+1自旋向上的值为+1自旋向下的值为-1。更进一步的简化是，

被认为是一个单位。图3显示了自旋最近邻居的相对位置，黑点被其周围邻居相互作用。图3，近邻耦合。位于(x,y)位置的黑点，正受到距离它一个晶格间距的最近点的相互作用。

为了最大化晶格边缘自旋的相互作用，它们被制成与晶格相对边缘的自旋相互作用。这被称为周期边界条件(periodic boundary condition)，如果我们考虑将二维格点折叠成一个自旋在该拓扑结构表面上的三维环面，就可以更好地表现出来。图4，一个三维环面的图解，它是具有周期边界条件的二维格点的代表。

1.4 计算问题

在一个辛模型的帮助下，我们可以继续为可观测的问题寻求解决方案。如果系统每一种可能状态的能量被指定，那么波尔兹曼分布函数，方程(1)，给出了系统处于每一种可能状态(在给定温度下)的概率，因此通过概率求和可以计算宏观量。这可以用磁化和能量作为例子来说明。对于任何固定状态，径向，当能量由哈密顿量(2)给出时，磁化强度与过量的向上或向下的自旋数成正比。

...(3)

M的期望值是

...(4)

E的期望值是

...(5)

这些计算从实际的角度提出了一个严重的问题。考虑到我们有两个自旋方向(向上和向下)，并且有N个自旋，这意味着有

个不同的状态。当N变大时，如果以这种方式计算，计算显然会成为一项令人生畏的任务。

使用计算机模拟来做这些计算似乎是一个很自然的建议，但是通过更仔细地研究方程(4)和(5)，很明显，使用这种方法在计算一个不可能的结果上浪费的计算精力与在一个很可能的结果上浪费的计算精力是一样多的。因此，一种更好的数字替代方法是使用模拟来生成关于“代表态”的数据。这些具有代表态构成了不同态的适当比例。这是有偏抽样的一种形式它本质上可以归结为满足以下条件

我们现在在一个更正式的环境中研究如何实现这一目标。

1.5 抽样和平均

可观测A(x)的热平均在正则系综中定义

...(6)

其中x是相空间中的一个向量，且

。配分函数Z由

而归一化的玻尔兹曼因子为

...(7)

这个概率给出了在热平衡中出现的结构x的实际统计权重。我们现在要考虑上述形式定义的离散情况。如果我们考虑相空间的有限部分，它将产生形式的平均值

...(8)

如果 在公式(8)中

，它会简化为式(6)。在相空间中取这种性质的简单样本的问题是，它不能保证概率分布在考虑的区域(不具有代表性)达到峰值。图5说明了这个问题。图5，这是一个产生以零为中心的高斯分布的简单采样示例，而关键数据在采样区域之外达到峰值。

因此，尝试一种更智能的采样技术来包括相空间中的重要区域是有意义的。我们想要一个选择点的过程，

，在相空间中有一个相关的概率

现在估计一个选定集合的热平均，

，将等式(8)简化为

...(9)

最明智的选择是

∝

。这个构造通过消去玻尔兹曼因子产生(9)的简单算术平均值，因此

...(10)

如果我们停下来思考我们在这一点上所要达到的目标，我们会发现我们正试图将无限相空间平衡时的概率分布减少到具有代表性的相空间有限集

的分布。问题是我们如何生成这个分布。

Metropolis等人通过引入连续状态

的马尔科夫过程的思想，改进了当时的传统思维，其中每个状态

都是由前一个状态

通过合适的转移概率

构造而成。为了成功地实现这一想法，必须建立一个详细的平衡方程，

...(11)

这个方程简单地说明，在平衡状态下，

和

的概率是相等的。如果我们现在取转移概率的比率，很明显，移动和反向移动只依赖于能量变化

。

...(12)

但是，这并不有助于惟一地指定W (xl xl)。为了做到这一点，我们引入

if

<0,

otherwise ... ...(13)

结果表明，利用该转移概率

，马尔可夫过程产生的状态分布

趋于均衡分布 在

。因此，当我们考虑相空间中更大的点集(

)时，该构造保持并以越来越高的精确度近似理论。

随时间变化的概率分布由马尔可夫主方程控制

...(14)

在

的热极限下，施加了详细的平衡方程(11)，结果如我们所料，

。此外，由于我们考虑的是一个有限系统，从逻辑上说，该系统是遍历的。

...(15)

如果我们将蒙特卡罗步数(mcs)视为时间单位，这种关系将以类似于前面讨论的算术平均值的方式减少。因此，我们面临的问题是，为了使平均时间与正则系综平均相等，系统是否遍历。如果我们把 mcs 看作时间的度量，那么这个条件就强加在系统上了。

...(16)

因此，Metropolis采样可以解释为在相空间中，由系统主方程(14)控制的沿随机轨迹的时间平均。

1.6 蒙特卡洛方法

事实是，我们能够成功地对抽样过程进行随机解释，这在允许我们引入蒙特卡罗方法作为一种技术来解决伊辛模型的可观测性方面起到了关键作用。

蒙特卡罗计算从根本上被定义为明确地使用随机变量，而随机过程具有随时间发展的随机过程的特征。因此，蒙特卡罗方法非常自然地适用于模拟发生随机过程的系统。根据我们所建立的关于采样类似于相空间中沿随机轨迹的时间平均的情况，可以使用蒙特卡罗方法来模拟这一过程。所采用的算法是围绕 Metropolis采样原则设计的。

从实现的角度来看，随机数的问题是这个过程的核心，它的成功取决于生成的数字是真正随机的这一事实。数值算法[9]广泛地处理了这个问题，并使用了本文的随机数生成器Ran1.c进行了模拟。这种与随机性的联系也是模拟这个名字的起源，因为蒙特卡洛的迷人的位置是运气和机会的同义词。

1.7 观测量的计算

感兴趣的观测量是

，

，

，

，

。

这些是按以下方法计算的

...(17)

采用上面的公式也可得

，

。

为了计算能量，我们使用方程(2)中给出的哈密顿量，

...(18)

引入1 / 2因子是为了考虑到旋转被计算了两次。式(18)同样用于确定

。

在居里温度下，我们预计这些量会有显著的涨落。一个很好的例子就是方差

。这使我们得到计算热容C和磁化率

的结论。

...(19)

...(20)

还计算了累积量。这将最终用来确定居里温度。

...(21)

c