熵的数值模拟方法与实操：态密度、热力学积分和纠缠熵

熵这个东西，在理论物理的语境里被讨论得太多，被真正算出来的太少。课本上写 (S = k_B \ln \Omega)，一行公式，干净利落；真到要拿数值说话的时候，绝大多数人第一步就卡住——你没法像测能量那样，在模拟里给熵加一个"观测量"。我这些年做统计物理和分子模拟相关的工作，凡是碰到自由能、相变、构象转变这类问题，最后都会绕回到同一件事上：熵的模拟。这篇就把这条链路从头到尾拆一遍，讲清楚熵在数值上到底怎么被"逼"出来，哪条路线适合什么场景，参数怎么定，坑在哪。

适读人群我大致分三类：刚接触蒙特卡洛、想算个比热和熵的入门者；做自由能计算、需要把熵从自由能里分离出来的从业者；以及做量子多体、需要计算纠缠熵的研究者。三类人的工具不同，但底层困境完全一致——熵是状态数的对数，而状态数在高维空间里是指数级的东西。理解了这一点，后面所有方法都只是"绕着指数做文章"的不同姿势。

1. 熵的模拟到底难在哪：一个绕不开的结构性困境

1.1 熵没有对应的力学量，只能通过"导数"或"计数"反推

能量、磁化强度、密度这些量，在模拟里都能写成位形坐标的瞬时函数，跑一步采样就能记一个数，最后取平均。熵不行。它不是一个关于微观状态的函数，而是一个关于系综的量，本质上描述的是"有多少个微观状态对应同一个宏观态"。你手上只有一条采样轨迹，轨迹里每个帧都是独立的微观状态，没有任何一帧"携带"熵的信息。

这是熵的模拟最反直觉的地方：它不是"算不出来"，而是"没有直接读取口"。所有可行的方法都只能走两条绕路。

第一条是走导数再积分。热力学里有一堆恒等式把熵对温度的导数跟可测量挂钩，比如 (C_v = T(\partial S/\partial T)_V)。既然比热可以靠能量涨落直接采样出来，那就能沿着温度把 (C_v/T) 积回去，倒推出熵。这条路线上的一切方法——热力学积分、自由能微扰、动力学方法——本质上都是"找一个熵的导数，测它，再积分"。

第二条是直接计数状态。既然熵是状态数的对数，那就想办法数。常规采样只在玻尔兹曼分布高概率的区域活动，低概率区域根本去不到，状态数数不全。于是就有了各种"强制平坦采样"的技巧，让随机游走在能量空间里均匀铺开，走一趟就把每一层的状态数记下来。

这两条路线的成本结构完全不同：积分路线精度高但需要做很多个温度点的独立模拟，计数路线一次跑完但收敛慢、参数敏感。选哪条，取决于你要的是"某个温度下熵的精确值"还是"整个温度区间的熵曲线"。

提示：如果你只需要两个状态之间的熵差（比如结合过程的构象熵损失），不要去做全温区积分，直接用自由能微扰或者热力学积分的窄区间版本，能把计算量压掉一个量级。

1.2 直接积分的维度灾难有多夸张

有人会问：状态数不就是把所有可能的微观状态列一遍数一数吗？我用一个具体的数来说明这条路为什么走不通。

拿最经典的二维自旋格子模型举例，取 (32\times32) 的格子，每个格点自旋取 ±1，总共有 (2^{1024}) 个状态。这个数大约是 (10^{308})。假设你有一台每秒列举 (10^{12}) 个状态的机器，从宇宙诞生算到今天，也只能覆盖其中微不足道的一小角。而这已经是最"温和"的模型了——换成连续坐标的分子体系，每个粒子的位置和动量都是连续变量，状态数是不可数无穷。

所以"数状态"这件事必须理解成对数密度：按能量分层，统计每一层能量有多少个状态，也就是态密度 (g(E))。能量是个连续（或准连续）量，分层之后层数大约是 (O(N)) 的量级，而不是指数级。这就把不可能变成了可能——前提是你有办法把每一层都采样得足够均匀。

注意：态密度不是"状态个数的近似"，它是精确的量，只是数值上有误差。分层宽度会引入离散化误差，这个误差不是随机的，而是系统性的，做高精度对比时要单独评估。

1.3 三类方法的适用范围，别拿错工具

我把常见的几条路线整理成一张表，这是选型的起点。很多人在项目里用错了方法，不是实现有问题，是路线本身不匹配需求。

选型上有个很实用的判定口径：如果体系有明显的相变或强涨落区，优先用能覆盖全温区的态密度方法；如果是液体或软物质，优先走分子动力学的两相热力学积分；如果是晶体低温性质，准谐振近似往往已经够用，不值得上重武器。

还有一个容易被忽视的点：熵的模拟结果几乎不可能独立验证。能量有精确解可对照的模型很少，绝大多数真实体系没有参照。因此实际项目里必须做交叉验证——用两种完全不同的方法算同一个体系的熵，误差在预期范围内才算过关。这个习惯我强烈建议从第一天就建立。

2. 先立一个能对照的靶子：二维自旋格子模型

2.1 模型定义与无量纲化，这一步省不得

做熵的模拟，我建议第一个练手项目永远选带精确解的模型。二维方格子自旋模型（每个格点一个 ±1 自旋，只跟最近邻耦合）是最合适的靶子：哈密顿量简单，有相变，有精确解，还有一大堆已知的临界指数可供对照。

哈密顿量写成 [ H = -J \sum_{\langle i,j\rangle} s_i s_j, \qquad s_i \in {+1, -1} ] 其中 (\langle i,j\rangle) 表示最近邻对。取周期性边界条件，避免边界自旋贡献不对称。

无量纲化是新手最容易偷懒、后果最严重的一步。我建议一开始就令 (J = 1)、(k_B = 1)，把所有量都用这两个单位表示。这样一来：

• 温度直接写成数值，临界温度理论值为 (T_c = 2/\ln(1+\sqrt{2}) \approx 2.269185)
• 能量、比热、熵全部以 (J) 和 (k_B) 为单位，输出时不用再乘系数
• 只有最后写论文或做实验对照时才乘回真实单位

我见过太多脚本里混着 (k_B) 和 (J) 的因子，跑出来的曲线形状对、数值差个 1.732，排查半天发现是某个地方多乘了一个 (\sqrt{3})。干脆不要给它们出现的机会。

提示：把单位常量定义在脚本顶部单独一行，比如J = 1.0; KB = 1.0，所有公式里只出现这两个变量名，永远不写裸数字。这一条看着啰嗦，能省掉后面大量的返工。

2.2 熵的两种来源：态密度与温度积分

在这个模型里，熵有两个等价但数值上完全不同的求法，正好用来互相印证。

第一条是态密度路线。把能量离散成若干个 bin，统计每个能量 bin 里的状态数 (g(E))，那么给定温度下的配分函数就是 [ Z(T) = \sum_E g(E), e^{-E/T} ] 有了配分函数，熵、自由能、内能、比热全部可以直接算出来： [ F = -T \ln Z,\quad U = \langle E\rangle,\quad S = (U - F)/T ] 这条路线的好处是一次采样拿到全部温度的结果，坏处是 (g(E)) 跨越几十个数量级，数值上很容易溢出，必须取对数存储、用对数域求和。

第二条是温度积分路线。从一个已知熵的参考温度出发，比如 (T \to \infty)，此时所有自旋完全随机，每个格点独立取两个值，总熵严格等于 (N\ln 2)（(N) 是格点数）。然后沿着温度往下降，每降一点就积一段： [ S(T) = S(\infty) + \int_T^{\infty} \frac{C_v(T')}{T'},dT' ] 这里 (C_v(T) = \beta^2 \mathrm{Var}(E))，(\beta = 1/T)，(\mathrm{Var}(E)) 是总能量的方差。用 (\beta) 作积分变量更稳，因为高 β（低温）区域的被积函数变化剧烈，而 β 网格可以做得很均匀。

两条路线算出来的熵曲线必须重合，这是最有价值的自检。

2.3 有限尺寸效应：小格子给出的熵到底是偏大还是偏小

有限尺寸这件事，做熵的模拟必须心里有数，因为它影响的方向跟能量、磁化不一样。

对 (L\times L) 的格子，总熵在 (T\to\infty) 时是严格的 (N\ln 2)，这个没误差。但在低温区，有限体系永远到不了真正的有序态——热涨落会把整个格子的磁化翻来翻去，形成两个简并的基态附近的小能隙。结果是低温熵不会趋于 0，而是趋于一个量级为 (\ln 2 / N) 的残余值。

对 (L=16)、(N=256) 的格子，这个残余大概是 (0.0027, k_B) 每自旋，看起来不大，但如果你做的是低温性质分析，这个残余会被误读成"基态简并度不为 1"。我一般在报告里明确写出"低温端的熵残余来自有限尺寸，非物理"，避免读者误判。

另一个更隐蔽的问题是临界区。临界点附近关联长度发散，(L) 必须远大于关联长度才能看到真实行为。判断标准是：如果你把 (L) 从 16 加到 32、再拉到 64，同一温度下的熵值变化超过统计误差，说明尺寸效应还没收敛，数据不能用。这个检查我每次都会做，比任何理论判据都直接。

3. 蒙特卡洛路线：怎么把比热积成一条熵曲线

3.1 单自旋翻转采样的实操细节

先从最基础的开始。单自旋翻转的接受判据是 [ \Delta E = 2 s_i \sum_{j \in \mathrm{nb}(i)} s_j, \qquad P_{\text{accept}} = \min\left(1, e^{-\Delta E / T}\right) ] 这是标准的 Metropolis 判据，写法上没有任何争议，但实现细节里埋着几个坑。

第一，每次扫描的翻转次数。一次"扫描"定义成 (L^2) 次翻转尝试，这样平均每个格点被尝试一次。不要用固定次数，否则不同尺寸的数据没法横向比较。

第二，预热的判据不能只看能量曲线。用能量自相关时间的两三倍作为预热长度是常见做法，但更稳的做法是同时监控能量和磁化强度绝对值，两者都稳定了才开始收录数据。临界点附近预热时间会拉长到几十倍，这是正常的，不要为了跑得快把预热砍掉。

第三，采样间隔与自相关。相邻扫描之间的能量是有相关性的，直接拿来算方差会低估 (C_v)。我的做法是每隔若干个扫描存一帧，间隔取自相关时间的 1 到 2 倍；或者更省事——存下全部数据，后面用分块法（binning）统一处理。分块法比硬性设定采样间隔更可靠，因为它直接告诉你统计误差。

注意：比热是从涨落里算的，涨落对相关性特别敏感。用未去相关的数据算 (C_v)，误差可以大到 30% 以上，而且这个误差是系统性的偏低，不是随机的。这是新手最常见的翻车点。

3.2 热力学积分：把 (C_v/T) 积回去的参数计算

具体的积分方案我是这么定的。以 (\beta) 为积分变量，(\beta) 从 0 积到目标值，公式为 [ S(\beta) = N\ln 2 - \int_0^{\beta} \frac{C_v(\beta')}{\beta'}, d\beta' ] 注意 (S(\beta=0) = N\ln 2)，积分号前是减号，因为 β 增大对应温度降低、熵减少。

被积函数 (C_v/\beta) 在 (\beta \to 0) 时的行为要特别处理：高温下 (C_v \propto \beta^2)，所以 (C_v/\beta \propto \beta \to 0)。这意味着 (\beta = 0) 这一点的被积函数取 0，而不是取 (C_v(0)/0)。如果不做这个处理，第一段的数值会炸掉。

网格密度怎么定？我的经验是：在临界区（(\beta \approx 0.44) 附近）用 0.005 的步长，其他区域用 0.02。梯形积分在这个密度下的截断误差大约在千分之几的量级，和统计误差相当。如果全部用均匀的 0.02 步长，临界峰附近的积分会引入大约 2% 的系统偏差，这个量级在我的验证里是能看出来的。

步长是不是够，有个很直接的检验方法：把网格加密一倍再算一遍，如果两次结果在统计误差内一致，说明截断误差可以忽略。这个检查做一次就够，不用每次都做。

3.3 可直接抄的代码骨架

下面算能量的函数，我用了周期性边界加np.roll的写法，比手写下标取模快不少，可读性也好。

import numpy as np def sample_energy(L, T, n_burn, n_keep, seed=20240501): """返回 n_keep 个总能量样本（J = 1, k_B = 1）""" rng = np.random.default_rng(seed) s = rng.choice([-1, 1], size=(L, L)).astype(np.float64) beta = 1.0 / T Es = np.empty(n_keep) for step in range(n_burn + n_keep): for _ in range(L * L): i = rng.integers(0, L) j = rng.integers(0, L) nb = (s[(i + 1) % L, j] + s[(i - 1) % L, j] + s[i, (j + 1) % L] + s[i, (j - 1) % L]) dE = 2.0 * s[i, j] * nb if dE <= 0.0 or rng.random() < np.exp(-beta * dE): s[i, j] = -s[i, j] if step >= n_burn: E = -0.5 * np.sum( s * (np.roll(s, -1, 0) + np.roll(s, 1, 0) + np.roll(s, -1, 1) + np.roll(s, 1, 1))) Es[step - n_burn] = E return Es

这里有个细节值得单独说：能量求和最后乘了 0.5，是因为每个键被数了两次。如果忘记这个 0.5，得到的能量是真实值的两倍，而且因为比热是二阶量，误差会被放大到四倍，但曲线形状完全正常，非常难发现。

积分部分：

def entropy_curve(L, betas, n_burn, n_keep): N = L * L integrand = np.zeros(len(betas)) for k, b in enumerate(betas): if b < 1e-9: continue E = sample_energy(L, 1.0 / b, n_burn, n_keep) Cv = b * b * E.var(ddof=1) # k_B = 1 integrand[k] = Cv / b # 梯形积分：从 beta=0 积到每个 beta S = np.empty(len(betas)) acc = 0.0 for k in range(1, len(betas)): acc += 0.5 * (integrand[k] + integrand[k - 1]) * (betas[k] - betas[k - 1]) S[k] = N * np.log(2.0) - acc S[0] = N * np.log(2.0) return S / N # 每自旋熵

这段代码返回的是每自旋熵，方便和理论值直接比。自检点很明确：函数值必须随 β 增大单调下降，且在 β=0 处严格等于 (\ln 2 \approx 0.693147)，低温端趋近于 0（有限尺寸下有几百个 (10^{-3}) 量级的残余）。

3.4 临界慢化：单自旋翻转什么时候必须换掉

单自旋翻转算法在临界点附近会严重退化。动力学临界指数在二维情况下大约是 2.17，意味着自相关时间随尺寸按 (L^{2.17}) 增长。(L=32) 时自相关时间可能是几十次扫描，(L=64) 就变成几百次，(L=128) 直接上四位数。这时候单靠加长采样时间已经撑不住了。

解决办法是换团簇算法。思路是把相邻的同向自旋按概率连成一片，整体翻转。连接概率取 [ p = 1 - e^{-2\beta} ] 这个概率的物理含义是：只有当两个相邻自旋同向时才会被连起来，连接概率正好抵消掉了界面能的贡献，使得团簇翻转的接受率恒等于 1。这是它最漂亮的地方——不需要算接受率，永远不会被拒绝。

def wolff_step(s, beta, rng): L = s.shape[0] i0, j0 = rng.integers(0, L), rng.integers(0, L) old = s[i0, j0] p_add = 1.0 - np.exp(-2.0 * beta) stack = [(i0, j0)] members = [(i0, j0)] seen = {(i0, j0)} while stack: i, j = stack.pop() for di, dj in ((1, 0), (-1, 0), (0, 1), (0, -1)): a, b = (i + di) % L, (j + dj) % L if (a, b) not in seen and s[a, b] == old and rng.random() < p_add: seen.add((a, b)) stack.append((a, b)) members.append((a, b)) for i, j in members: s[i, j] = -s[i, j] return len(members)

团簇算法的动力学指数降到 0.25 左右，临界点的采样效率提升是数量级的。代价是它对某些观测量（比如需要逐点操作的量）不那么友好，但算能量、磁化这些全局量完全没问题。我的习惯是：远离临界区用单自旋翻转，靠近临界区切团簇算法，两种算法在同一温度点跑一遍对比，结果一致就放心用。

4. 直接数状态数：平坦直方图类方法的实操口径

4.1 态密度的核心思路与收敛判据

前面说过，状态数是按能量分层的密度。平坦直方图类方法的核心动作就是：人为修改采样概率，让随机游走在能量空间里走成一条近似水平的直线，而不是被玻尔兹曼因子压成指数衰减。这样每个能量层都能被访问到足够多次，密度就能估出来。

具体做法是维护一张能量层的访问直方图 (H(E)) 和一张态密度的对数表 (\ln g(E))。每次尝试翻转，先算新能量，如果这一层被访问过的次数多，就压低接受概率。迭代一段时间后，直方图趋于平坦，就把它清零、把修改幅度减半，进入下一轮。修改幅度从 1 开始，按 0.5、0.25 这样逐步减小，直到小到可以忽略。

收敛判据我一般用两条：一是修改幅度降到 (10^{-8}) 以下，二是最后几轮的直方图相对起伏小于 5%。只满足第一条不够，因为直方图起伏大说明还没真正铺平，只是"看起来"收敛了。

有个很实际的问题：这个方法天然是串行的（每一步依赖上一步的密度表），并行化很难。小体系还行，(32\times32) 大概跑几十亿次翻转，普通工作站要跑上一两天。大体系基本没法用，这时候还是老老实实回去做温度积分。

4.2 参数选择：分 bin 宽度和修改幅度

分 bin 的宽度决定了离散化误差。太宽，态密度曲线的细节被抹平，相变位置会偏移；太窄，很多 bin 长时间无人访问，直方图起伏大，收敛慢。

经验参数是：总能量范围除以 200 到 400，作为 bin 宽度。二维 (32\times32) 格子，基态能量大约 -2048，最高能量大约 +2048，跨度 4096，分成 300 个 bin，每个 bin 宽约 13.6。这个宽度下，临界区的精细结构能分辨出来，同时大部分 bin 都能被有效访问。

另一个参数是"什么时候把直方图清零"。常见做法是每一轮访问次数达到 (10^4) 到 (10^5) 次量级就清零。太早清零，密度表还没累积够信息；太晚清零，低访问率的 bin 拖慢整体进度。

提示：跑完之后一定要检查态密度表里有没有"空洞"——某些 bin 的访问次数几乎为零。如果有，通常是能量范围取错了，或者初始构型太特殊导致某些层根本去不到。空洞不处理，后面的配分函数求和会直接算错。

4.3 两种方法的交叉验证怎么做

拿到态密度之后，配分函数在每一个温度都能算出来，熵曲线一次成型。这时候最有价值的工作是做交叉验证。

具体做法：选三到五个温度点，比如 (T = 1.0, 1.8, 2.269, 3.0, 5.0)，在这些点上用态密度算一次熵，再用温度积分算一次熵。两者之差如果落在各自统计误差的一倍以内，说明两条路线都对。如果某个温度点偏差明显，通常说明这个温度附近态密度支撑不足，要回去检查那个能量段的分 bin 和采样。

我做过一次统计，在 (32\times32) 体系上，两种方法的平均偏差大概在 0.5% 量级，临界点附近稍大一些，能到 1%。这个量级可以作为参考，但不同实现肯定有差异，关键是建立"自己的实现能复现到什么程度"这个基准线。

还有一点：态密度方法得到的熵是绝对熵，直接就是物理量；温度积分也得到绝对熵，前提是高温参考点取对。两边都是绝对量，才能直接比。

5. 分子动力学的路子：从自由能里把熵剥出来

5.1 熵的分解方式决定了误差结构

到了分子体系，问题复杂度上一个台阶。这时候的主流做法不是直接算熵，而是先算自由能，再把内能从自由能里减掉：(S = (U - F)/T)。所以本质上还是"算自由能 + 算内能"。

对溶液、蛋白这类体系，常用的分解是把熵拆成几块：

这个分解很重要，因为它告诉你误差主要来自哪一块。绝大多数情况下，构象熵是主要不确定度来源，平动和振动项的贡献加起来往往还不如构象项的一半。所以如果项目时间紧，应该把计算资源压在构象熵上，而不是在振动分析上抠小数。

5.2 热力学积分的实操：耦合参数怎么设计

自由能的热力学积分，核心是引入一个耦合参数 (\lambda)，从 (\lambda=0)（参考态，没有相互作用）连续变到 (\lambda=1)（全相互作用），然后积 [ \Delta F = \int_0^1 \left\langle \frac{\partial H}{\partial \lambda} \right\rangle_\lambda d\lambda ] 积分变量是 (\lambda)，被积函数是每个 (\lambda) 下的系综平均。

这里最关键的经验是：(\lambda) 网格不能均匀取，必须在变化剧烈的地方加密。(\partial H/\partial \lambda) 通常在 (\lambda \to 0) 附近出现尖峰，因为这时候一个粒子刚从"无相互作用"变成"有相互作用"，周围什么都没有，能量变化很突兀。均匀取 11 个点，低 λ 区域的积分误差能占到总分量的 10% 以上。

我的做法是分成三段处理：低 λ 段用软核势把相互作用平滑地"打开"，避免奇点；中间段均匀取点；高 λ 段适当加密。这样每段的被积函数都比较平，积分稳定。

还有个大坑：(\lambda) 相邻点之间的构型必须能互相转换。如果两个 λ 值对应的构型差太多，中间的采样会卡住，出现"滞后"现象。解决办法是让相邻 λ 的差足够小，或者从后一个 λ 的构型开始前一个 λ 的采样，用这种串行方式保证连续性。

5.3 粗粒化带来的"熵口径"问题

做粗粒化模型的时候，有个概念上的坑必须提前想清楚：粗粒化会把一部分自由度"积分掉"，这部分自由度对应的熵会被吸收到有效势里。所以粗粒化模型算出来的熵，和全原子模型算出来的熵，不是同一个东西。

举个具体例子：把四个原子的基团简化成一个珠子，那么珠子之间的有效相互作用已经包含了原来基团内部构象的平均效果。用这个有效势算出来的构象熵，天然比全原子模型小，差值恰好是被积分掉的那部分内部自由度的贡献。

这个差值不是误差，是定义差异。问题在于很多人拿粗粒化的熵和全原子的实验值直接比，发现对不上，然后开始怀疑代码。正确的做法是：要么两边都粗粒化，要么明确报告"这里的熵不含某某自由度的贡献"。我见过论文里因为这个问题被审稿人抓住反复修改的情况。

6. 量子体系：纠缠熵是另一套坐标系

6.1 约化密度矩阵，一步都不能错

量子体系里的熵，最常见的对象是纠缠熵：把系统分成 A、B 两部分，对 B 求迹得到 A 的约化密度矩阵 (\rho_A)，然后 [ S_A = -\mathrm{Tr}\left(\rho_A \ln \rho_A\right) = -\sum_i w_i \ln w_i ] (w_i) 是 (\rho_A) 的本征值。形式上很简单，实现上最容易错的是子系统的划分和态矢的索引顺序必须严格对应。

一个长度为 (L) 的自旋链，态矢有 (2^L) 个分量。如果约定"第一个格点对应最高位"，那么前 (L_A) 个格点构成的子系统对应态矢最前面的 (L_A) 个比特位。这时候把态矢 reshape 成 ((2^{L_A}, 2^{L-L_A})) 的矩阵，直接把矩阵乘自己的共轭转置，就得到 (\rho_A)。

问题在于很多人写代码时用的约定和建哈密顿量时的约定不一致。行索引和列索引的方向反了，算出来的纠缠熵数值看起来完全正常（正数、量级合理），但物理意义全错。这个错误特别难发现，因为熵的数值不会"异常"，只是不对。

import numpy as np def entanglement_entropy(psi, L, LA): """约定：格点 0 对应态矢最高位""" psi = psi.reshape(2 ** LA, 2 ** (L - LA)) rho = psi @ psi.conj().T w = np.linalg.eigvalsh(rho) w = w[w > 1e-12] # 去掉数值负值和零 return float(-np.sum(w * np.log(w)))

这个函数里那个w > 1e-12的截断是必须的，因为数值求本征值会给出微小的负值，取对数直接变 NaN。

6.2 精确对角化的验证流程

希尔伯特空间维数是 (2^L)，能精确对角化的上限大概在 (L=16) 到 (L=20)（(2^{20}) 约一百万维，稠密矩阵需要几个 GB 内存）。这个尺寸只够做验证，做不出有意义的物理结论，但验证的价值极大。

我建议的验证步骤是这样的：

第一步，用一个解析可解的例子校准。最简单的选择是无相互作用的自旋链或者完全二聚化的极限，这两种情况下纠缠熵都有闭式结果。跑一遍，看数值和解析式差多少，这个差值就是整条链路的"本底精度"。

第二步，检查两个极限：子区为空时纠缠熵必须严格为 0；子区为整个系统时，如果总态是纯态，纠缠熵也必须为 0。这两个检查能抓住大部分索引错误。

第三步，检查对称性：如果没有破坏某种对称性，那么互补子系统的纠缠熵必须相等，即 (S_A = S_B)。这一条对纯态严格成立，是抓索引混乱的利器。

做完这三步，才有资格去算真正想看的东西。

6.3 面积律和数值陷阱

纠缠熵的标度行为跟系统尺寸有关：有能隙的基态，纠缠熵随子区边界长度增长（面积律），一维情况下趋于饱和；临界体系的基态，一维情况下纠缠熵随子区尺寸对数增长。这个差别在数值上是很清楚的信号，可以用来判断体系处于什么相。

数值陷阱主要有两个。

一是基态简并。如果基态有简并，精确对角化给出的本征矢可能是简并子空间里任意一个叠加态，这个态本身可能没有物理意义，算出来的纠缠熵会随实现细节乱跳。处理办法是：先看本征值谱，确认第一激发态的能量间隙，如果间隙小于数值精度，就要用对称性把简并子空间分开，或者干脆加上一个小的对称破缺项。

二是子区尺寸的奇偶性。半填充附近，奇数和偶数子区尺寸的结果可能系统性地不同，这在做标度拟合时会造成假象。我的做法是至少取两种奇偶性的尺寸，画在同一张图上，看曲线是否衔接自然。

7. 常见问题与排查实录

7.1 问题速查表

这些年带人做熵的模拟，遇到的问题高度集中，我整理成一张表，对号入座能省不少时间。

7.2 几个不容易被发现的小坑

第一个是随机数质量。熵的模拟对随机数敏感度远高于能量模拟，因为熵的量级是几十到几百，误差是累积的。用低质量的线性同余发生器，跑出来的熵曲线会有肉眼可见的抖动，而且换种子后结果会变。我的标准做法是固定用同一个高质量生成器、固定种子、跑三条不同种子的轨迹，三条结果在误差内一致才采信。如果换种子结果就变，那不是统计误差，是系统性问题。

第二个是单位制的最后一步。前面用无量纲单位跑完全部计算，最后乘回真实单位的时候，熵要乘 (k_B)、能量要乘 (J)。如果 (J) 本身是通过实验数据拟合出来的，那么拟合的不确定度会直接传递到熵上。这个传播关系很多人不算，直接在论文里报一个精确到小数点后三位的熵值，实际上有效数字可能只有一位。

提示：报熵的结果时，永远同时说明"以什么为单位"和"误差来源构成"。绝对熵的不确定度往往比熵差大一个量级，因为绝对熵要把高温参考点到目标点的积分全走一遍，误差全部累积。

第三个是积分路径的选取。热力学积分的路径不唯一，换一条路径得到的结果必须相同（状态函数）。这既是自检手段，也是优化手段：有些路径上的被积函数起伏小、积分好算。我一般在正式跑之前，先用一两条廉价路径比一比，选起伏小的那条走。

第四个是采样的"陷阱"。有些体系存在能量亚稳态，单温度模拟会被困在里面出不来，表现为熵曲线在某温度处出现台阶。这时候必须做多温度并行（不同温度之间交换构型），或者用平坦采样把整个能量空间铺开。判断方法很简单：从不同初始构型出发跑同一温度，结果不一致，就是有陷阱。

8. 我个人踩过的几条经验

做熵的模拟，最大的体会是：精度不是靠算得多，是靠自检多。我见过太多脚本，代码写得漂亮、跑得飞快，但从来没做过高温极限自检，结果曲线整体偏移一个常数，最后被审稿人一句"高温极限都不对"打回来。熵这个量天生没有直观图像，错了不会"看起来怪"，只能靠一系列硬性的检查点兜住。我现在的习惯是每个熵计算脚本里至少写五个断言：高温熵、低温行为、单调性、双方法一致性、不同种子一致性。这五条过了，才敢往报告里放数。

另一个体会是关于误差棒。熵的误差不像能量那么好估，因为它是积分量，误差会沿着积分路径累积。我的处理办法是把每个温度点的误差单独算出来，然后做误差传播积分，同时再用"整条路径重跑三次取散布"的办法交叉验证。两种估计差太多，说明某个环节有系统性偏差。这个过程比较繁琐，但它决定了你报的数是"计算结果"还是"可信结果"。

最后一个实际建议：如果你刚开始做这个方向，不要一上来就搞大体系或者真实材料。找一个有精确解的模型，把整条链路走通，把五个断言全部写出来，把两种方法都对上，把误差结构摸清楚。这个练手过程大概需要一到两周，之后换任何体系都是把模块替换掉的事。反过来，如果跳过这一步直接上真实体系，你连结果对不对都不知道，那才是真正的浪费。