由于本人水平不足，本文难免会有错误，如发现错误或有任何不懂的欢迎在评论区指出

本部分主要参考：

很多分子模拟的教材和其他的计算物理教材都有蒙特卡洛模拟方法的介绍。

通过前两部分（0-1，0-2）的介绍，我们只需要得到体系状态的概率分布：

就可以通过统计的方法求得物理量。但是其中配分函数𝑄𝑁𝑉𝑇并不容易计算，并且学习过统计物理会知道，如果我们知道了配分函数的解析表达式，其实可以直接通过配分函数来计算相关物理量而不需要使用（0-2）中的公式。

这里使用蒙特卡洛模拟的方法生成符合此概率分布的体系状态，从而避免了直接计算配分函数，这样就需要使用（0-2）中公式来计算物理量了。

蒙特卡洛（Monte Carlo，MC）模拟原理

一般来说，蒙特卡洛模拟是按照一定概率，随机地将体系从状态 𝑆𝑖 转换到状态 𝑆𝑗。转换概率为 𝑊𝑖𝑗，且要求 𝑊𝑖𝑗 只与当前状态 𝑆𝑖 和下一个状态状态 𝑆𝑗 有关（与之前的过程无关，即转换是无记忆的），即所谓马尔科夫过程（Markov process）。

一系列马尔可夫过程构成马尔可夫链（Markov chain），在平衡时有关系（全局平衡条件，global balance equations）：

其中 𝑃(𝑆𝑖) 为体系在 𝑆𝑖 时的概率。

引用 Thijssen 的比喻，将每个状态 𝑆𝑖 视为一个个水桶，而 𝑃(𝑆𝑖) 为每个水桶中的水量，𝑊𝑖𝑗 𝑃(𝑆𝑖) 为从水桶 𝑖 流向水桶 𝑗 的水量，从而平衡时应有（对每一个 𝑖 ）：

流向水桶 𝑖 的总水量（左边）=从水桶 𝑖 流出的总水量（右边）

这里模拟时仅需满足充分条件即可，即细致平衡条件（detailed balance condition）：

继续上述的比喻，则这里认为每两个水桶之间相互流动的水量是相等的（自然满足上述平衡条件）。

相比全局的平衡条件，细致平衡条件没有考虑成环的概率流动等其他形式。上图的例子不满足细致平衡条件。

这里我们需要构造 NVT 系综下的超额部分的体系状态分布，即要求（详见0-1的“NVT 系综下的进一步化简”，将物理量分为理想部分和超额部分）：

代入上式得到：

选取满足上述的 𝑊𝑖𝑗 构造马尔科夫链，体系概率稳定时即为所需生成的概率分布，由此避免了直接计算配分函数。

为了方便构造，将 𝑊𝑖𝑗 拆分成 𝑊𝑖𝑗 = 𝑃𝑖𝑗𝐴𝑖𝑗，令其中 𝑃𝑖𝑗 为选择概率，𝐴𝑖𝑗 为接受概率，故有：

对应模拟中，选择概率 𝑃𝑖𝑗 意为在状态 𝑆𝑖 中选择到状态 𝑆𝑗 的概率，而 𝐴𝑖𝑗 为选择了上述状态后接受这个选择的概率。注意在模拟中即使没有接受也应计入统计。

Metropolis 算法

这里选取对称的选择概率：

而接受概率选择：

（代入满足上式）

即当下一个状态（𝑆𝑗）的能量小于当前状态（𝑆𝑖）时一定接受这个选择，而当下一个状态的能量大于于当前状态时则有 exp(-βΔE) 的概率接受这个选择。

这种选取选择概率和接受概率的算法即为 Metropolis 算法

补充

由于蒙特卡洛方法只是一个统称，许多只要涉及随机选择的算法都可以称为蒙特卡洛方法。故对于本文中使用蒙特卡洛方法生成马尔科夫链得到概率分布用于统计的方法，一般称为马尔科夫链蒙特卡洛算法（Markov Chain Monte Carlo, MCMC），或者 Metropolis Monte Carlo 算法。

当然这里主要用于模拟物理体系，所以笔者称其为蒙特卡洛模拟，后续介绍动力学蒙特卡洛模拟（Kinetic Monte Carlo，KMC）时也称这种方法为原始的蒙特卡洛模拟。

总结

本部分介绍了蒙特卡洛模拟的原理和方法，在下一部分（1-2）将会用其模拟单原子气体（跳过 Ising 模型部分）。

如发现错误或有任何不懂的欢迎在评论区指出