计算统计物理：蒙特卡罗方法与分子动力学

§32.1为什么需要计算机模拟

对相互作用系统，配分函数一般无法解析计算（第17、18章）。直接数值求和也不可能：$10\times10$ 的伊辛模型就有 $2^{100}\approx10^{30}$ 个组态；$N$ 个粒子的位形积分是 $3N$ 维的积分，用每维只取 10 个点的网格也需要 $10^{3N}$ 个点。解决办法是随机抽样：不去遍历所有组态，而是按照玻尔兹曼分布抽取有代表性的组态，用样本平均代替系综平均。这有两种方式：

• 蒙特卡罗方法：构造一个随机过程，使它的平稳分布就是玻尔兹曼分布（梅特罗波利斯、罗森布鲁斯夫妇、特勒夫妇，1953）；
• 分子动力学：直接对牛顿方程作数值积分，按各态历经假设（§2.7）用时间平均代替系综平均（奥尔德与温赖特，1957）。

§32.2重要性抽样

要计算 $\langle A\rangle = \sum_iA_ie^{-\beta E_i}/Z$。若均匀地随机选取组态，绝大多数被选中的组态能量很高、玻尔兹曼因子小得可以忽略，对平均几乎没有贡献——高维空间中，玻尔兹曼分布集中在相空间里一个极小的区域（§2.3.4、§4.4）。重要性抽样：直接按概率 $P_i = e^{-\beta E_i}/Z$ 抽取 $M$ 个组态 $i_1,\dots,i_M$，则

若样本相互独立，由中心极限定理（提示 C6），$\bar A$ 的统计误差为 $\sqrt{\mathrm{Var}(A)/M}$，按 $M^{-1/2}$ 减小，与系统的维数无关。困难在于：$Z$ 未知，无法直接按 $P_i$ 抽样。

§32.3梅特罗波利斯算法

思路：构造一个马尔可夫过程（§23.6），使它的跃迁速率满足细致平衡 (23.18)：

由 §23.6，这样的过程从任何初态出发，分布都会趋向玻尔兹曼分布（只要过程是各态历经的：任何组态都能经过有限步到达任何其他组态）。注意 (32.2) 只涉及能量差，不需要知道 $Z$。

梅特罗波利斯规则：

1. 从当前组态 $i$ 出发，提出一个试探组态 $j$（例如翻转一个随机选取的自旋，或把一个随机选取的粒子移动一小段随机距离），提议是对称的：从 $i$ 提出 $j$ 的概率等于从 $j$ 提出 $i$ 的概率；
2. 计算 $\Delta E = E_j - E_i$；
3. 以概率 $\min\left(1,e^{-\beta\Delta E}\right)$ 接受 $j$，否则保留 $i$（并把 $i$ 再计一次）。

验证细致平衡：设 $E_j>E_i$，则 $W(i\to j)\propto e^{-\beta(E_j - E_i)}$，$W(j\to i)\propto1$（提议概率相同，相消），两者之比正是 (32.2)；$E_j<E_i$ 时同理。

能量降低的步总被接受，能量升高的步以玻尔兹曼因子的概率被接受：系统主要在低能组态附近活动，又能借助热涨落越过势垒。

§32.4例：二维伊辛模型

对伊辛模型，翻转自旋 $s_i$ 的能量变化只依赖于它的近邻：$\Delta E = 2Js_i\sum_{j\in\text{近邻}}s_j$，在正方晶格上只能取 $0,\pm4J,\pm8J$。把格点像棋盘一样分成两套，同一套中的自旋互不相邻，可以同时更新。下面是一个完整的 Python 程序（以 $J = k_{\mathrm B} = 1$ 为单位，周期性边界）：

import numpy as np

def ising_mc(L, T, n_sweep=6000, n_therm=1500, seed=1):
    rng = np.random.default_rng(seed)
    s = np.ones((L, L), dtype=np.int8)               # 从全部向上出发
    i, j = np.indices((L, L))
    sublattices = [(i + j) % 2 == c for c in (0, 1)] # 棋盘的两套格点
    m_samples = []
    for sweep in range(n_sweep):
        for mask in sublattices:
            nb = (np.roll(s, 1, 0) + np.roll(s, -1, 0)
                  + np.roll(s, 1, 1) + np.roll(s, -1, 1))
            dE = 2 * s * nb                          # 翻转每个自旋的能量变化
            accept = rng.random((L, L)) < np.exp(-np.maximum(dE, 0) / T)
            s[mask & accept] *= -1                   # 梅特罗波利斯规则
        if sweep >= n_therm:                         # 先弃去趋向平衡的部分
            m_samples.append(abs(s.mean()))
    return np.mean(m_samples)

print(ising_mc(32, 2.0))   # 约 0.911

用这个程序得到的每个自旋平均磁化的绝对值 $\langle\lvert m\rvert\rangle$，与杨振宁的严格结果 (§18.6) 比较（$k_{\mathrm B}T_{\mathrm c}/J = 2.269$）：

远离 $T_{\mathrm c}$ 时，即使很小的晶格也给出精确的结果；在 $T_{\mathrm c}$ 附近与 $T_{\mathrm c}$ 以上，有限尺寸效应很明显：有限的系统中没有真正的相变（§21.4），$\langle\lvert m\rvert\rangle$ 随 $L$ 增大才缓慢地趋于零。§32.5 将说明如何利用这种尺寸依赖。

§32.5统计误差、临界慢化与有限尺寸标度

关联时间。马尔可夫链中相继的样本不是独立的：每一步只改变一个或几个自旋。若某个量的自关联函数（§22.4，以"扫描"即每个自旋平均被尝试一次为时间单位）的积分关联时间为 $\tau$，则 $M$ 次扫描中大约只有 $M/(2\tau)$ 个独立样本，统计误差相应地增大。

临界慢化。在临界点附近，关联长度 $\xi$ 很大，局域的单自旋翻转要改变大小为 $\xi$ 的区域需要很多步：$\tau\propto\xi^z$，在 $T_{\mathrm c}$ 处 $\tau\propto L^z$。对二维伊辛模型的梅特罗波利斯算法，$z\approx2.17$。斯文森–王（1987）与沃尔夫（1989）的集团算法一次翻转一整团相互关联的自旋（集团的选取规则保证细致平衡），把 $z$ 降低到 0.3 以下，使临界点附近的模拟效率提高了许多个数量级。

有限尺寸标度。有限系统中 $\xi$ 不能超过 $L$。由 $\chi\propto\lvert t\rvert^{-\gamma}\propto\xi^{\gamma/\nu}$，在尺寸为 $L$ 的系统中可以写成标度形式

（当 $L\gg\xi$ 时 $\Phi(x)\propto\lvert x\rvert^{-\gamma}$，回到无限系统的结果；当 $\xi\gg L$ 时 $\chi$ 被截断在 $L^{\gamma/\nu}$。）由此：$\chi$ 的峰高 $\propto L^{\gamma/\nu}$（二维伊辛模型为 $L^{7/4}$）；峰的位置偏离 $T_{\mathrm c}$ 的距离 $\propto L^{-1/\nu}$；把不同 $L$ 的数据画成 $\chi L^{-\gamma/\nu}$ 对 $L^{1/\nu}t$，应落在同一条曲线上。宾德累积量 $U_L = 1 - \frac{\langle m^4\rangle}{3\langle m^2\rangle^2}$ 是一个无量纲的量，在 $T_{\mathrm c}$ 处与 $L$ 无关，所以不同 $L$ 的 $U_L(T)$ 曲线交于一点，给出 $T_{\mathrm c}$ 的精确估计。今天三维伊辛模型最精确的临界温度，以及精度仅次于共形自举（§20.5）的临界指数，大多来自这类有限尺寸标度分析。

§32.6分子动力学

韦尔莱算法。对每个粒子的牛顿方程 $m\ddot{\mathbf r} = \mathbf F(\mathbf r)$ 作数值积分。把 $\mathbf r(t\pm\Delta t)$ 泰勒展开（提示 A3）：$\mathbf r(t\pm\Delta t) = \mathbf r\pm\mathbf v\Delta t + \frac12\mathbf a\Delta t^2\pm\frac16\dddot{\mathbf r}\Delta t^3 + O(\Delta t^4)$，两式相加，奇次项相消：

（韦尔莱，1967。）等价的"速度韦尔莱"形式为：$\mathbf r(t + \Delta t) = \mathbf r + \mathbf v\Delta t + \frac{\mathbf F(t)}{2m}\Delta t^2$，$\mathbf v(t + \Delta t) = \mathbf v + \frac{\mathbf F(t) + \mathbf F(t + \Delta t)}{2m}\Delta t$。这个算法是时间可逆的（把 $\Delta t$ 换成 $-\Delta t$ 就沿原路返回，与 §P1.2 一致），并且是一种"辛"积分方法：像真实的哈密顿演化一样保持相空间的几何结构（其中包括相空间体积，刘维尔定理，§2.6）。正因为如此，总能量的误差在长时间内只是有界地振荡而不会系统地漂移，这对长时间模拟至关重要。时间步长通常取最快振动周期的约 1/20：对伦纳德–琼斯氩约 10 fs，对含 C–H 键的分子约 1 fs。

测量（周期性边界条件下的 $N$ 个粒子）：

• 温度：由均分定理，$\frac32Nk_{\mathrm B}T = \langle K\rangle$（$K$ 为总动能）；
• 压强：由 (17.16) 的推导推广到一般的力，$p = nk_{\mathrm B}T + \frac{1}{3V}\left\langle\sum_{i<j}\mathbf r_{ij}\cdot\mathbf F_{ij}\right\rangle$；
• 对关联函数 $g(r)$：统计粒子对的距离直方图（§17.7），再傅里叶变换得到 $S(k)$（§22.3），可以与衍射实验直接比较；
• 扩散系数：由均方位移 $\langle\lvert\mathbf r(t) - \mathbf r(0)\rvert^2\rangle\to6Dt$，或由格林–久保公式 (23.9)；黏度、热导率同样可以用格林–久保公式（§24.7）由平衡模拟算出。

恒温。上述算法保持能量守恒，抽样的是微正则系综。要模拟恒定温度，可以给每个粒子加上朗之万方程中的摩擦与随机力，两者的比例由 (23.6) 确定，这样就抽样正则系综（朗之万恒温器）；也可以用确定性的诺泽–胡佛方法。

一个历史上的发现。奥尔德与温赖特（1957）以及伍德与雅各布森（1957）的模拟发现：只有排斥、没有任何吸引的硬球，在足够高的密度下也会结晶。有序的晶体反而比无序的流体熵更高——在高密度下，排列整齐的球各自有更大的"活动空间"。这个纯粹由熵驱动的相变在当时出人意料，后来在胶体悬浮液中被实验证实。

§32.7自由能的计算

自由能与熵不是某个量的系综平均，不能直接从样本平均得到，需要专门的方法。

热力学积分。设哈密顿量依赖于参数 $\lambda$（例如从理想气体 $\lambda = 0$ 逐渐"打开"相互作用到 $\lambda = 1$）。由 $F = -k_{\mathrm B}T\ln Z$，$\frac{\partial F}{\partial\lambda} = -\frac{k_{\mathrm B}T}{Z}\frac{\partial Z}{\partial\lambda} = \frac1Z\sum\frac{\partial\mathcal H}{\partial\lambda}e^{-\beta\mathcal H}$，即

在若干个 $\lambda$ 值上做平衡模拟，测量 $\langle\partial\mathcal H/\partial\lambda\rangle_\lambda$，再数值积分。

维多姆粒子插入法。化学势 $\mu = F(N+1) - F(N) = -k_{\mathrm B}T\ln(Z_{N+1}/Z_N)$。由 (17.2)，$\frac{Z_{N+1}}{Z_N} = \frac{Q_{N+1}}{(N+1)\lambda^3Q_N}$；而 $Q_{N+1}/Q_N = V\langle e^{-\beta\Delta U}\rangle_N$，其中 $\Delta U$ 是在 $N$ 粒子系统中随机位置插入一个"试验粒子"时它与其他粒子的相互作用能，平均对 $N$ 粒子的平衡组态与插入位置进行。于是

第一项就是 (3.17)。对稠密液体，随机插入的粒子几乎总与别的粒子重叠，$e^{-\beta\Delta U}\approx0$，这个方法就失效了，需要更复杂的技巧。远离平衡的拉伸实验与模拟中，还可以用雅津斯基等式（§33.4）由非平衡功的分布求出平衡自由能差。

§32.8本章小结

1. 高维的系综平均必须用重要性抽样；统计误差按样本数的 $-1/2$ 次方减小。
2. 梅特罗波利斯算法用满足细致平衡的马尔可夫链抽取玻尔兹曼分布，只需要能量差。
3. 临界点附近有临界慢化 $\tau\propto L^z$，集团算法可以克服它；有限尺寸标度把尺寸效应变成测量临界指数的工具。
4. 分子动力学用时间可逆、保持相空间体积的韦尔莱算法积分牛顿方程；输运系数可以由平衡模拟的格林–久保公式求出。
5. 自由能需要专门的方法：热力学积分、粒子插入法、非平衡功的等式。

自测题

1. 证明：若梅特罗波利斯算法的提议不对称（从 $i$ 提出 $j$ 的概率 $q_{ij}\ne q_{ji}$），把接受概率改为 $\min\left(1,\frac{q_{ji}}{q_{ij}}e^{-\beta\Delta E}\right)$ 仍满足细致平衡（梅特罗波利斯–黑斯廷斯算法）。
2. 修改 §32.4 的程序，计算每个自旋的磁化率 $\chi = \beta N(\langle m^2\rangle - \langle\lvert m\rvert\rangle^2)$，对 $L = 16$、$32$ 找出峰的位置与高度，并检验峰高之比是否接近 $2^{7/4}\approx3.4$。
3. 用速度韦尔莱算法积分一维谐振子（$\omega = 1$），取 $\Delta t = 0.1$，验证能量误差不随时间累积；再用简单的欧拉法 $(x,v)\leftarrow(x + v\Delta t,\ v - x\Delta t)$（同时更新，右边都用旧值）比较。［答：欧拉法的能量每步乘以 $1 + \Delta t^2$，不断增长；若先更新 $x$、再用新的 $x$ 更新 $v$，就成为"辛欧拉法"，能量误差有界］
4. 由 (32.5) 证明：若 $\mathcal H(\lambda)$ 是 $\lambda$ 的线性函数，则 $F(\lambda)$ 是 $\lambda$ 的凹函数。［提示：计算 $\partial^2F/\partial\lambda^2 = -\beta\left(\langle(\partial_\lambda\mathcal H)^2\rangle - \langle\partial_\lambda\mathcal H\rangle^2\right)\le0$］