在模拟零场Ising模型的时候，我们给定一系列温度，通过Monte Carlo方法在每个温度下达到平衡，求得序参量平均磁矩$m_0$的值，从而给出平均磁矩随温度变化的曲线，从中观察$m_0$发生突变（从零到一个有限值）的温度，最终得到这一相变温度。这是一种直接模拟的思路。

然而直接模拟法需要十分注意滞后效应。Ising模型的模拟中，簇算法能够很好地帮助系统翻越自由能垒，从而克服了滞后效应。但是簇算法有着比较强的专一性，所以在某些体系下并不容易找到这样的算法。另一种思路是并行退火，其目的也是翻越自由能垒。并行退火有着一定的普适性。

原理

热力学积分一种研究固液相变的常用方法。相变点处，两相共存并且系统处于平衡态。对于微正则系综（孤立体系，NVE），由热力学第二定理知熵$S$最大时系统处于平衡态；对于正则系综（恒温恒容封闭体系，NVT），Helmholtz自由能$F$最小时系统处于平衡态；对于恒温恒压系综（NPT），Gibbs自由能$G$最小时系统处于平衡态。通常Monte Carlo模拟正则系综，这里以正则系综为例。

假如我们希望知道在某个温度下两相（密度不同，给定）是否能够稳定共存（相平衡），那么只需要分别计算其自由能。只要自由能相等，我们认为这是一个相平衡点（相变点）。

但是自由能（以及熵）并不能直接通过模拟得到。因为与相空间的全体有关，而模拟不可能遍历全体相空间。

事实上这并不成问题，因为我们只要设定一个参考系统，比较自由能变$\Delta F$即可。由热力学第二定律可知，我们可以在目标点和参考点之间找一个可逆过程，通过热力学积分来求自由能变。

首先，参考系统是一个自由能已知的系统，例如理想气体和低温简谐晶体。在计算液体的自由能时通常取理想气体作为参考系统。

其次，所谓可逆过程，只要在积分路径上系统随时达到平衡并且热力学函数连续即可，也就是说不能跨一级相变。而且，在模拟中，选择非物理的某条路径也可以，比如人为设定一个参数：

$\begin{equation} U=(1-\lambda) U_1+\lambda U_2 \end{equation}$

这里$U$为过程中系统的势能，随着参数$\lambda$从0连续变化到1，$U$从参考系统的势能$U_1$变化为目标系统的势能$U_2$。在正则系综下，自由能变化量为

$\begin{equation} F(\lambda=1)-F(\lambda=0)=\int _ {\lambda=0}^{\lambda=1} d\lambda\Bigl<\frac{\pd U(\lambda)}{\pd \lambda}\Bigr>_\lambda \end{equation}$

若是跨一级相变，如气液相变，则考虑从临界点（二级相变点）外侧绕过去，如图所示

对于液相，计算非常直接，液体的Helmholtz自由能$F$为

$\begin{equation} \frac{F(\rho)}{Nk_BT}=\frac{F^\text{id}(\rho)}{Nk_BT}+\frac{1}{k_BT}\int_0^\rho d\rho'\left(\frac{P(\rho')-\rho k_BT}{\rho'^2}\right) \end{equation}$

积分可以分两步：第一步，由高温理想气体（临界温度以上）等温可逆压缩到目的密度；第二步，等密度（等容）可逆降温到目的温度。沿着积分路径求数值积分即可。其中$P(\rho’)$的计算肯定需要使用模拟的方法（MC或MD）。

对于固相，从理想气体开始积分是不可能的，因为不存在一个“临界点”绕过去。但是固体的自由能可以从一个已知自由能的理想固体模型求出。所以，重点是设计是个恰当的、自由能可精确求解析解或数值解的理想固体模型，

参考资料

Frenkel, Daan, and Berend Smit. Understanding molecular simulation: from algorithms to applications. Vol. 1. Elsevier, 2001.
Thermodynamic integration. Wikipedia (2018.5) （援引 J Kästner; et al. 文章，称：Umbrella sampling is a related free energy method. It adds a bias to the potential energy. In the limit of an infinite strong bias it is equivalent to thermodynamic integration.）
Zhang, Kai, and Patrick Charbonneau. “Monte Carlo approach for studying microphases applied to the axial next-nearest-neighbor Ising and the Ising-Coulomb models.” Physical Review B 83.21 (2011): 214303.