计算物理中的熵模拟：态密度、配分函数与量子熵可复现流程

1. 从研究场景切入：熵的模拟到底在算什么

熵的模拟，在理论物理研究里不是把一条公式抄进代码那么简单。它更像是在给一个系统的“可能性空间”做人口普查：哪些微观构型存在、各占多少权重、温度变化时这些权重如何重新分配。做得好，你能从数值上看到相变附近自由能如何弯曲、比热如何起峰、纠缠如何增长；做得糙，画出来的曲线也能“看起来像那么回事”，但熵的常数项、能量格点、采样误差会一起埋进结果里，最后连自己都不敢拿它去解释物理。这篇复盘围绕“熵的模拟”展开，重点讲清楚三件事：热力学熵如何从态密度里算出来，量子熵如何从密度矩阵或张量网络里取出来，以及一套能落地复现的采样流程该怎么搭。适合正在做计算物理练习的研究生、需要给理论模型补数值证据的从业者，也适合刚接触蒙特卡洛、想把“配分函数”从公式变成可操作对象的高年级本科生。下文会以二维伊辛模型和一个小型量子自旋系统为例，但方法本身可以迁移到更一般的统计模型与量子多体问题。

1.1 为什么理论物理研究绕不开熵

在理论物理里，能量往往是最先被盯上的量，因为哈密顿量写出来就有；熵则麻烦得多，它不是一个局域算符的期望值，而是整个相空间或希尔伯特空间里概率分布的全局性质。也正因为麻烦，熵的模拟才有价值。微正则系综里，熵直接对应态密度：(S(E)=k_B\ln g(E))，其中 (g(E)) 是能量为 (E) 的微观构型数。正则系综里，熵又和配分函数绑在一起：(F=-k_BT\ln Z)，(S=(U-F)/T)。你只要能把 (g(E)) 或 (Z) 算准，自由能、比热、熵变、相变强度就有了共同的落脚点。反过来，如果只算平均能量，很多细节会被掩盖：两个系统可以有相近的 (U(T))，但熵和自由能差别很大，临界行为也可能完全不同。

我在实际项目里最深的体会是，熵的模拟最考验“定义是否闭合”。你得先回答：我要的是热力学熵、冯·诺依曼熵、Rényi 熵，还是轨迹熵产生？这几种熵长得像，数值路线却差很远。热力学熵通常从态密度或配分函数走，量子熵通常从密度矩阵走，纠缠熵还要先划分子系统。定义没闭合就开跑，后面所有参数选择都会摇摆。比如有人拿经典 Wang-Landau 采样去套量子模型，却忘了量子模型没有简单的“能量构型计数”，这时候 (g(E)) 的物理含义已经变了。先把熵的定义、系综和模型自由度对齐，是后面所有步骤的第一道门槛。

另一个容易被低估的点是参考态。离散自旋模型的态密度可以绝对计数，因为总构型数是 (2^N)，高温极限给出一个天然归一化条件：(\sum_E g(E)=2^N)。连续模型或粗粒化模型就没这么幸运，熵里常带一个与分辨率有关的常数。这个常数不影响比热、也不影响相对熵，但会影响你报告“绝对熵”时的可信度。所以我在写模拟方案时，会先明确哪些量是平移不变的，哪些量必须靠外部参考点固定。把这件事提前说清楚，比事后解释“曲线整体偏移是正常的”要省很多麻烦。

1.2 熵的模拟到底在模拟什么

把问题拆到底层，熵的模拟其实在估计一个权重分布。经典模型里，这个分布是 (p(E)\propto g(E)e^{-\beta E})；量子模型里，是密度矩阵 (\rho=e^{-\beta H}/Z) 的本征值分布。经典模拟的难点在于 (g(E)) 跨越几十个数量级，直接均匀采样能量空间会严重偏向高能区；量子模拟的难点在于希尔伯特空间随粒子数指数增长，精确对角化很快撞墙。两条路线的共同点是：不能靠“多跑几步”硬扛，必须设计有偏采样或压缩表示，把有限算力花在真正重要的区域。

以二维伊辛模型为例，(N=L^2) 个自旋，能量范围从 (-2JN) 到 (+2JN)，单自旋翻转改变能量为 (4J) 的整数倍。这个离散结构非常适合做态密度采样。可如果你用普通 Metropolis 在低温跑，系统会被困在低能谷里，高能区几乎访问不到；在高温跑，低能区又采样不足。Wang-Landau 的思路是反过来：不按玻尔兹曼权重走，而是让随机游走在能量轴上尽量平坦，用直方图平坦度作为收敛判据，逐步修正态密度。这个“让能量空间被均匀访问”的设计，是经典熵模拟里最核心的一步。

量子侧则更像线性代数加统计推断。小系统可以直接构造哈密顿量矩阵，做精确对角化，得到本征值和本征矢，再算热态密度矩阵或基态约化密度矩阵。系统一大，就得换量子蒙特卡洛、张量网络或神经网络量子态。这里的关键不是“方法名字够不够新”，而是你能不能控制误差来源：蒙卡有符号问题，张量网络有键维截断，神经网络有优化不收敛。熵的模拟在量子侧尤其怕“看起来收敛”：能量可能已经稳了，但熵对低概率本征值很敏感，截断一改，数值就漂。所以我在量子熵项目里会同时盯三个量：能量、配分函数、以及熵本身随截断参数的变化趋势，三者不一致就不急着下结论。

1.3 哪些人适合看这篇复盘

如果你正在做统计物理模型的数值练习，想从“会跑 Metropolis”进阶到“会算自由能和熵”，这篇内容会对口。如果你在做量子多体的小系统精确对角化，想确认密度矩阵熵的数值细节，也能直接拿走代码骨架。如果你只是听说过“熵的模拟”这个词，想弄明白它和普通分子动力学、普通蒙卡有什么区别，前面的场景拆解和后面的方法对比表也足够建立框架。需要的前置知识不多：会一点 Python、知道配分函数是什么、能看懂矩阵本征值分解即可；更深的理论推导我不会铺开，重点放在参数怎么定、误差怎么查、结果怎么交叉验证。

我也得提前泼一盆冷水：熵的模拟没有“一键出图”的通用脚本。不同模型的能量格点、边界条件、对称性、序参量都不一样，参数换个模型就得重调。下文给的参数是基于二维伊辛模型和小型量子自旋系统的常见实践，属于可复现的起点，不是放之四海皆准的固定配方。真正做研究时，你必须用自己的模型做收敛性测试，至少用两种独立方法交叉验证。比如经典侧用 Wang-Landau 算 (g(E))，再用普通 Metropolis 算 (U(T)) 对一遍；量子侧用精确对角化算小系统，再和高温展开或张量网络在可重叠区域比对。这个交叉验证的习惯，比任何单个参数都重要。

2. 方案选型：不同熵定义对应不同模拟路线

2.1 热力学熵路线：态密度与配分函数是主线

经典统计模型里，最稳妥的熵模拟路线是先拿态密度 (g(E))，再由它生成整条温度曲线。原因很直接：(g(E)) 一旦算准，配分函数 (Z(\beta)=\sum_E g(E)e^{-\beta E}) 就在一次求和里完成，平均能量 (U=\sum_E E,g(E)e^{-\beta E}/Z)、比热 (C_v=(\langle E^2\rangle-\langle E\rangle^2)/(k_BT^2))、熵 (S=k_B(\ln Z+\beta U)) 都能顺带得到。你不需要对每个温度单独跑一轮模拟，这对研究相变特别友好，因为相变附近温度点很密，逐点跑蒙卡既慢又难保证误差一致。态密度路线的代价是采样难度高：能量空间的低端和高端访问概率差得离谱，必须用平坦直方图类方法或并行回火来补。

我选 Wang-Landau 作为主线，不是因为它最新，而是因为它的输出形态最适合讲熵。它直接给 (\ln g(E))，误差传播路径清楚，收敛判据也能可视化。相比之下，普通 Metropolis 加多次模拟退火也能估自由能，但要把整条熵曲线拼出来，热力学积分会引入路径依赖和积分误差。热力学积分不是不能用，它在分子模拟里非常成熟：(S(T)=S(T_0)+\int_{T_0}^{T} C_v/T,dT)。可对离散自旋模型，既然态密度可以离散计数，就没必要绕这个弯。选型的逻辑是：模型能量空间离散且维度不高，优先态密度；模型连续、自由度多，才考虑热力学积分或自由能微扰。

注意：态密度方法的前提是能量空间能被有限个格点覆盖。连续模型直接套 Wang-Landau 会碰到格点宽度选择问题，格点太粗会抹平精细结构，太细又会让直方图平坦化变得极慢。连续系统更常见的是先做离散化近似，或者改用并行回火加直方图重加权。

2.2 量子熵路线：密度矩阵、纠缠熵与张量网络

量子侧的熵模拟要先分清热熵和纠缠熵。热熵对应热态 (\rho=e^{-\beta H}/Z) 的冯·诺依曼熵 (S=-{\rm Tr}\rho\ln\rho)，它描述的是有限温度下混合态的不可知程度。纠缠熵则通常指基态或某个本征态下，把系统分成 A、B 两部分后约化密度矩阵 (\rho_A={\rm Tr}_B|\psi\rangle\langle\psi|) 的熵。两者数值路线有交集，但物理含义不同：热熵随温度升高而增大，纠缠熵在量子相变附近才可能出现显著增长，且强烈依赖分区方式。做模拟前不把这两个概念分开，后面很容易把“温度升高导致熵增”和“纠缠增长导致熵增”混在一张图里。

小系统最可靠的办法是精确对角化。粒子数 (N) 不大时，希尔伯特空间维度是 (2^N) 或 (3^N)，构造稀疏矩阵，用 Lanczos 求低能本征对，热态密度矩阵可以由本征分解直接拼出来。这个方法的优点是误差可控，缺点是维度墙来得很快：(N=16) 时维度是 65536，稠密矩阵已经吃不消，必须用稀疏存储；(N=20) 以上就要非常小心。系统再大，张量网络是更现实的选择：矩阵乘积态（MPS）适合一维或准一维，PEPS 适合二维但代价高，量子蒙特卡洛适合大尺度但受符号问题限制。我的选型习惯是：先把小系统精确解当“标尺”，再用张量网络或蒙卡去够大尺度，最后在两者重叠的参数区间对熵做交叉验证。

量子熵模拟还有一个隐蔽的坑：密度矩阵的本征值动态范围极大。热态低温下，低能本征值占主导，高能本征值小到 (10^{-15}) 以下，直接算 (\ln\rho) 会出现 (0\cdot\infty) 和负本征值的数值噪声。我的处理方式是先做本征分解，把本征值截断到正数，用 (x\ln x) 在 (x\to0) 时趋于 0 的极限单独处理小本征值。这个细节看起来琐碎，但它决定低温熵曲线是否平滑。张量网络侧则要盯键维截断：键维从 32 加到 64，熵如果还在明显漂移，就说明截断没收敛，不能拿当前结果去拟合临界指数。

2.3 方法对比：别拿锤子找钉子

这张表不是排名，而是选型检查清单。我见过最常见的错误是：模型明明是二维阻挫量子磁体，却直接上量子蒙特卡洛，跑了两天发现符号问题严重，样本有效数少得可怜；或者模型只是 16 个自旋的小系统，却非要搭张量网络，调键维和规范的时间比精确对角化还长。选型时先问三个问题：能量空间是否离散、希尔伯特空间维度是否可控、目标熵是热熵还是纠缠熵。答案清楚了，方法基本就定了。剩下的才是参数调优，而不是一开始就在方法名字上纠结。

3. 核心细节解析：从态密度到熵的可复现流程

3.1 能量离散化与能量格点设计

经典自旋模型的能量格点不是随便分的。以二维方晶格伊辛模型为例，哈密顿量 (H=-J\sum_{\langle ij\rangle}s_i s_j)，(s_i=\pm1)，总键数是 (2L^2)，所以能量范围是 ([-2JN,+2JN])，其中 (N=L^2)。单自旋翻转时，能量改变 (\Delta E=2J s_i\sum_{j\in{\rm nbr}(i)}s_j)。在方晶格上，邻居和只能取 (-4,-2,0,2,4)，乘以 (2J s_i) 后，(\Delta E) 是 (4J) 的整数倍。这意味着能量格点步长取 (4J) 就够了，再细就是浪费直方图 bin。这个推导看着简单，但它决定了 bin 数和采样效率。步长取错，要么漏掉可达能量，要么把多个能量塞进一个 bin，后者会让熵曲线出现虚假平台。

具体算一下：(L=16)，(N=256)，(E_{\min}=-512J)，(E_{\max}=512J)，步长 (4J)，bin 数是 ((512-(-512))/4+1=257)。(L=32) 时，(N=1024)，bin 数是 1025。这个规模对内存完全不是问题，真正的成本在采样步数。每个 bin 最终阶段至少要有足够多的访问次数，平坦化才不会靠偶然涨落过关。我的经验是：最终修改因子降到 (\ln f\sim10^{-7}) 时，每个 bin 的有效访问量最好在 (10^4) 到 (10^5) 量级；总步数随系统变大而增加，(L=16) 用 (10^7) 到 (10^8) 步比较常见，(L=32) 往往要更多，且要结合多独立 run 平均。这个数字不是定理，但能作为一个起跑线。

注意：边界条件会改变能量格点结构。周期性边界和开放边界的键数不同，缺陷或无序模型可能让能量步长不再整齐。每次换模型，先写一个小脚本枚举所有可达能量或做短程随机游走，把实际能量集合打印出来，再决定 bin 映射。别直接沿用伊辛模型的步长。

3.2 Wang-Landau 采样参数怎么定

Wang-Landau 的核心是维护态密度 (g(E)) 和访问直方图 (H(E))。随机游走在能量空间进行，从旧能量 (E_{\rm old}) 提议新能量 (E_{\rm new})，接受概率取 (\min(1,g(E_{\rm old})/g(E_{\rm new})))，接受后更新 (g(E_{\rm new})\leftarrow g(E_{\rm new})f)，并让 (H(E_{\rm new})) 加一。当直方图足够平坦，把修改因子开方，重置直方图，进入下一轮。这个流程里，初始修改因子 (f)、平坦度阈值、终止条件、能量窗是四个关键参数。初始 (f) 常取 (e) 或 2，取 (e) 意味着第一轮对态密度的修正比较温和；平坦度阈值常取 0.8，意思是所有已访问 bin 的直方图计数不低于平均值的 80%。阈值太高会让每轮拖得很久，太低则容易过早进入下一轮，留下系统偏差。

终止条件通常看 (\ln f)，降到 (10^{-6}) 到 (10^{-8}) 之间停。停得太早，(g(E)) 的低概率区域噪声大；停得太晚，收益递减。我的做法是至少跑两条独立链，用不同随机种子，比较 (\ln g(E)) 的逐 bin 差异。如果两条链在高能端和低能端的差异都落在可接受范围，再合并；如果只有中间能量一致、两端打架，说明能量窗边缘采样不足，需要扩展窗或分段采样。能量窗本身也要先探：用短程 Metropolis 从随机态和高低温态各跑一段，记录能量范围，再把窗稍微放宽一点。窗太窄会截断真实分布，窗太宽会让大量 bin 永远访问不到，平坦度判据直接失效。

这里还有一个实现细节：接受概率用比值 (g(E_{\rm old})/g(E_{\rm new}))，不要用 (g(E_{\rm new})/g(E_{\rm old}))。我见过不止一次因为写反导致态密度越跑越偏，最后高温端熵出现负值。另一个细节是更新顺序：先算接受概率，接受后再更新 (g) 和 (H)。如果把更新放在接受判断之前，相当于偷偷改变了提议分布，收敛判据会失真。代码短不代表可以省检查，写完先用 (L=4) 的小系统暴力枚举所有构型，把精确 (g(E)) 和模拟结果对一遍，这个基准测试半小时就能做完，后面能省几天排查。

3.3 从 (\ln g(E)) 到热力学熵的换算与参考点

Wang-Landau 输出的是相对态密度，整体乘一个常数不确定。对离散伊辛模型，可以用高温极限定标：(Z(\beta=0)=\sum_E g(E)=2^N)。于是把模拟得到的 (g_{\rm sim}(E)) 乘以 (2^N/\sum_E g_{\rm sim}(E))，就得到绝对计数的态密度。这个定标步骤很关键，因为微正则熵 (S(E)=k_B\ln g(E)) 里的常数项会直接影响绝对熵数值。若只关心比热、相对熵或自由能差，不定标也能用，但报告绝对熵时必须有参考点。连续模型没有 (2^N) 这种天然计数，通常要选一个参考温度，用热力学积分或已知极限固定常数。

有了 (g(E))，正则系综的量可以一次性算出来。取 (k_B=1)，配分函数 (Z(\beta)=\sum_E g(E)e^{-\beta E})，平均能量 (U=\sum_E E g(E)e^{-\beta E}/Z)，比热 (C_v=(\langle E^2\rangle-\langle E\rangle^2)/T^2)，熵 (S=\ln Z+\beta U)。计算时先对能量做平移，把 (e^{-\beta E}) 里的最大指数提出来，避免低温下溢。比如令 (E_0=\min E)，计算 (e^{-\beta(E-E_0)})，最后再把 (E_0) 的贡献补回去。这个技巧不改变物理，但能救回低温曲线。我通常还会检查两个极限：(T\to\infty) 时 (S/N\to\ln2)；(T\to0) 时 (S) 应趋于基态简并熵。两个极限对不上，先别怀疑物理，先查定标和数值下溢。

3.4 误差、自相关与有限尺寸标度

熵的误差不像平均能量那样直观，因为 (g(E)) 的误差会通过对数放大。我的做法是保留多条独立 run，对 (\ln g(E)) 逐 bin 做平均和标准差，再用 bootstrap 重采样生成 (S(T)) 的置信带。不要只跑一条长链然后把样本切成几段当独立样本，Wang-Landau 的修改因子在早期是非平稳的，切段估计会偏乐观。最终阶段固定 (g(E)) 后，可以再做一轮普通 Metropolis 采样来估计热力学量的自相关时间，用 (N_{\rm eff}=N/(2\tau+1)) 折算有效样本数。这个折算对 (C_v) 特别重要，因为比热是涨落量，自相关影响比平均值大。

有限尺寸效应也要提前规划。二维伊辛模型的临界温度精确值是 (T_c=2J/(k_B\ln(1+\sqrt2))\approx2.269J/k_B)。你跑 (L=8,12,16,24)，比热峰位会随 (L) 移动，峰高也会变。想从熵曲线提取临界行为，至少要做三到四个尺寸，并用有限尺寸标度拟合。别拿单个 (L=16) 的曲线去断言临界指数。我的习惯是先看熵密度 (S/N) 随尺寸是否收敛，再看比热峰高是否按预期增长，最后才做标度拟合。如果 (S/N) 在小尺寸和大尺寸之间明显漂移，先查边界条件和能量窗，别急着上拟合。

4. 实操过程：二维伊辛模型熵模拟完整演示

4.1 环境与数据结构

这套演示用 Python 加 NumPy 就够了，不需要复杂框架。核心数据结构只有四个：自旋矩阵spins、态密度对数ln_g、直方图hist、能量到 bin 的映射。自旋矩阵用(L,L)的整型数组，取 (\pm1)；能量用整数或浮点都可以，但索引时必须统一成整数 bin。为了避免反复计算总能量，我维护当前能量E，每次翻转只算局部能量差。这个增量更新是性能关键，(L=32) 时全量重算能量会让循环慢一个数量级。随机数生成器固定种子，方便复现；但正式统计要多 seed 跑独立链，不能只靠一个种子。

import numpy as np L = 16 N = L * L J = 1.0 rng = np.random.default_rng(20240517) spins = rng.choice([-1, 1], size=(L, L)) E_min = -2 * J * N E_max = 2 * J * N dE = 4 * J n_bins = int((E_max - E_min) / dE) + 1 def total_energy(s): e = 0.0 for i in range(L): for j in range(L): e -= J * s[i, j] * ( s[(i+1) % L, j] + s[(i-1) % L, j] + s[i, (j+1) % L] + s[i, (j-1) % L] ) return 0.5 * e def energy_index(e): return int(round((e - E_min) / dE)) E = total_energy(spins) ln_g = np.zeros(n_bins, dtype=float) hist = np.zeros(n_bins, dtype=np.int64)

这段代码里total_energy除以 2 是因为每条键被两个自旋各数了一次。energy_index用round而不是直接取整，是为了防止浮点误差把边界能量分错 bin。实际跑大系统时，我会把邻居表提前建好，把二维索引展平成一维，减少 Python 循环开销。演示代码保留二维写法是为了可读性，正式跑之前建议先做向量化或改用 Numba，否则 (L=16) 以上会明显变慢。数据结构不复杂，但每个索引映射都要用暴力枚举验证一遍，尤其是能量边界。

4.2 核心更新与接受准则

Wang-Landau 的更新循环可以写得很短，但每个分支都要有明确意图。随机选一个自旋，算翻转前后的能量差，映射到 bin，计算接受概率。接受就用新能量更新当前状态，然后修正态密度和直方图；不接受就保持原状态，但当前 bin 的直方图是否计数要一致。常见做法是无论接受与否，都对当前能量 bin 计数一次，这样直方图反映的是访问频率，不是接受频率。平坦度判据只看已访问 bin，未访问 bin 不参与比较。修改因子降到阈值后停止，输出 (\ln g) 和直方图。

flatness = 0.8 ln_f = 1.0 min_ln_f = 1e-7 while ln_f > min_ln_f: hist[:] = 0 steps = 0 while True: i = rng.integers(0, L) j = rng.integers(0, L) s = spins[i, j] nb_sum = ( spins[(i+1) % L, j] + spins[(i-1) % L, j] + spins[i, (j+1) % L] + spins[i, (j-1) % L] ) dE_local = 2.0 * J * s * nb_sum E_new = E + dE_local old_idx = energy_index(E) new_idx = energy_index(E_new) # 接受概率用旧态密度比新态密度 accept_prob = min(1.0, np.exp(ln_g[old_idx] - ln_g[new_idx])) if rng.random() < accept_prob: spins[i, j] = -s E = E_new idx = new_idx else: idx = old_idx ln_g[idx] += ln_f hist[idx] += 1 steps += 1 visited = hist > 0 if steps > 1000 and visited.sum() > 0: h_mean = hist[visited].mean() if hist[visited].min() >= flatness * h_mean: break ln_f *= 0.5

这段循环里，ln_g存的是对数态密度，更新时直接加ln_f，避免大数相乘溢出。接受概率用exp(ln_g[old] - ln_g[new])，如果指数太大就截断到 1。平坦度判断放在steps > 1000之后，是为了防止刚开始直方图太少导致误判。实际跑的时候，我会把每轮steps、ln_f、已访问 bin 数打印出来，观察是否出现“平坦度很快达标但只访问了少数 bin”的假收敛。如果能量窗内存在永远访问不到的区域，直方图会一直不包含那些 bin，平坦度判据反而容易通过，这是最隐蔽的坑之一。

注意：ln_f每轮减半是常见做法，但不必死守。有人用 (f\leftarrow\sqrt{f}) 的对数形式，本质一样。关键是最终ln_f要足够小，且每轮直方图重置后要重新跑够步数。不要在一轮还没平的时候就手动降ln_f，那相当于人为制造偏差。

4.3 运行记录与结果检查

一次练习性质的复现里，我用 (L=8,12,16) 分别跑了两条独立链，初始 (f=e)，平坦度 0.8，终止 (\ln f=10^{-7})。(L=16) 的单链在普通笔记本上大约几分钟到十几分钟量级，具体取决于实现是否向量化。跑完后先看三张图：(\ln g(E)) 是否关于 (E=0) 近似对称、直方图最终是否覆盖整个能量窗、以及由 (g(E)) 算出的 (U(T)) 是否和普通 Metropolis 结果重合。对称性来自伊辛模型的自旋翻转对称 (E\to-E)，如果模拟结果明显不对称，先查能量窗和 bin 映射，再查随机数。直方图覆盖不全，说明能量窗或采样步数有问题。

结果检查里最有用的一步是交叉验证平均能量。用 Wang-Landau 得到的 (g(E)) 算 (U(T))，再用普通 Metropolis 在几个温度点独立跑，比较两条曲线。两者在相变附近可能有小幅差异，但整体趋势必须一致。如果 (U(T)) 在低温端偏离，常见原因是低能区采样不足或 bin 太粗；如果高温端偏离，常见原因是能量窗上限不够或定标常数错了。比热曲线则看峰位是否接近 (T_c\approx2.269J/k_B)，峰高是否随 (L) 增大。这个检查不需要高深理论，但能快速暴露大多数实现错误。

我还习惯保存原始 (\ln g(E)) 和直方图，而不是只存最终曲线。因为后面想换温度点、换定标方式、做 bootstrap，都要回到原始数据。只存一张 (S(T)) 图，等于把中间证据扔了。多链结果合并时，不要简单平均 (\ln g)，先检查每条链的逐 bin 差异，再决定是平均、加权还是丢弃异常链。异常链往往在低能端或高能端出现大幅偏移，直接平均会把偏差带进最终结果。

4.4 量子小系统：密度矩阵熵的数值验证

量子侧先用小系统建立标尺。以两个自旋 1/2 的海森堡二聚体为例，构造哈密顿量矩阵，做本征分解，算有限温热态密度矩阵 (\rho=e^{-\beta H}/Z)，再算冯·诺依曼熵 (S=-{\rm Tr}\rho\ln\rho)。这个例子小到可以手算校验，适合检查代码路径。密度矩阵本征值可能出现数值上的微小负值，先截断到正数，再用 (x\ln x) 的极限处理零本征值。温度趋近零时，熵应趋近基态简并熵；温度趋近无穷时，熵应趋近 (\ln D)，(D) 是希尔伯特空间维度。两个极限都对上，才说明熵的计算没有漏项。

import numpy as np sx = np.array([[0, 1], [1, 0]], dtype=float) sy = np.array([[0, -1j], [1j, 0]], dtype=complex) sz = np.array([[1, 0], [0, -1]], dtype=float) I2 = np.eye(2) def kron(a, b): return np.kron(a, b) # 两个自旋的 Heisenberg 二聚体，J=1 H = (kron(sx, sx) + kron(sy, sy) + kron(sz, sz)) beta = 1.0 w, V = np.linalg.eigh(H) w0 = w.min() rho = V @ np.diag(np.exp(-beta * (w - w0))) @ V.conj().T rho /= np.trace(rho) evals = np.linalg.eigvalsh(rho) evals = np.clip(evals, 0.0, None) mask = evals > 1e-15 S_vN = -np.sum(evals[mask] * np.log(evals[mask])) print("eigenvalues:", evals) print("S_vN:", S_vN) print("S_max:", np.log(len(evals)))

这段代码的重点不是模型多复杂，而是数值路径清楚：本征分解、热权重、归一化、截断、熵求和。实际做纠缠熵时，先把基态写成矩阵形式，对 B 子系统做偏迹得到 (\rho_A)，再重复同样的熵计算。偏迹实现要特别小心索引顺序，二维系统分区时更容易错。我的做法是先用一个可分离态验证：如果 A、B 没有纠缠，(\rho_A) 应是纯态，熵接近 0；再用一个最大纠缠态验证，熵应接近 (\ln2)。这两个基准过了，再上真实模型。量子熵的代码往往不长，错都错在索引和截断上。

5. 常见问题与排查技巧实录

5.1 典型故障速查表

这张表里最值得展开的是“直方图平坦但只覆盖部分 bin”。这个现象很阴，因为平坦度判据只比较已访问 bin，未访问 bin 不参与，所以程序会告诉你“已经平坦”，实际上能量空间根本没走完。解决办法不是简单调高平坦度阈值，而是先确认能量窗内哪些能量可达。对离散模型，可以从小系统枚举推断；对连续模型，可以做短程随机游走统计。确认可达集合后，把不可达 bin 从判据里排除，或者把能量窗缩到可达范围。否则你会在一个假收敛结果上继续跑，越跑越自信。

另一个高频问题是“低温熵负值”。这通常不是物理，而是计算 (e^{-\beta E}) 时下溢成 0，导致 (Z) 变成 0 或极小，取对数后出问题。处理方式是在求和前减去最小能量 (E_0)，算 (Z'=\sum g(E)e^{-\beta(E-E_0)})，最后用 (\ln Z=\ln Z'-\beta E_0) 补回。这个技巧对所有配分函数计算都适用。类似地，量子密度矩阵先减基态能量再归一化，也能避免低温下整矩阵下溢。

5.2 踩坑心得：那些文档里不写的细节

第一个心得是：先做小系统暴力基准，再跑大系统。很多人一上来就 (L=32)，跑完发现曲线不对，回头查代码成本极高。(L=4) 的伊辛模型只有 16 个自旋，总构型 65536 个，暴力枚举几分钟就能得到精确 (g(E))。拿它对照 Wang-Landau 输出，接受准则、bin 映射、定标、熵换算全都能验证。这个习惯我在每个新模型上都会保留，哪怕模型看起来很简单。基准测试不是浪费时间，它把“物理错误”和“实现错误”分开了。

第二个心得是：能量窗宁可先窄后宽，不要一开始铺太大。窗太宽会让大量 bin 永远访问不到，平坦度判据失真，采样效率也被稀释。我的做法是先用短程 Metropolis 从随机态和高低温态各跑一段，记录能量分位数，再把窗设到覆盖 99.9% 访问范围，边缘留少量余量。跑完一轮后看直方图两端是否接近零，如果边缘 bin 计数很高，说明窗可能截断了分布，需要扩窗重跑。这个流程比拍脑袋设窗稳得多。

第三个心得是：不要用一条链的最终结果直接报误差。Wang-Landau 的误差估计没有普通 Metropolis 那么直接，多链 bootstrap 是最省事的方案。至少三条独立链，逐 bin 算均值和标准差，再重采样生成 (S(T)) 的置信带。如果三条链的 (\ln g) 在某个能量段分叉，说明那段采样不可靠，先别做物理结论。误差棒不是装饰，它决定了你能不能区分两个相近的临界指数。

第四个心得是：量子熵的截断参数要写成可扫描的配置，不要硬编码。张量网络的键维、精确对角化的本征值截断、量子蒙卡的 Trotter 步长，都会影响熵。把这些参数放进配置表，跑一组扫描，看熵随参数的变化是否收敛。我见过能量已经收敛到小数点后六位、熵还在漂移的情况，因为熵对低概率本征值更敏感。只盯能量会误判收敛。

5.3 结果可信度自检清单

在把熵曲线放进报告之前，我会过一遍自检清单。第一，高温极限对不对：经典伊辛模型 (S/N\to\ln2)，量子系统 (S\to\ln D)。第二，低温极限对不对：熵应趋近基态简并熵，而不是负值或突然跳水。第三，热力学一致性对不对：由 (S) 求导得到的 (C_v) 和直接涨落算的 (C_v) 是否一致，(F=U-TS) 是否和 (Z) 算的自由能一致。第四，尺寸趋势对不对：熵密度随 (L) 增大是否趋于稳定，比热峰是否按预期移动。第五，独立方法对不对：至少用一个更简单的方法在部分参数区交叉验证。

这份清单里，热力学一致性最容易被跳过。有人只算 (S(T))，不回头检查 (C_v=T(\partial S/\partial T))，结果熵曲线看着平滑，导数却出现负比热。负比热不一定错，但在稳定平衡态里出现就要警惕。我的做法是把 (S(T)) 拟合成平滑样条再求导，和涨落公式对比；两者在误差范围内一致才继续。这个检查能抓住定标错误、bin 太粗、温度网格太稀等问题。熵的模拟不怕算得慢，就怕错得安静。

6. 扩展方向与个人经验

6.1 从平衡态熵到动态熵产生

平衡态熵算稳之后，很容易想往非平衡走。一个自然的扩展是轨迹熵产生：对一条随机轨迹 (x(t))，比较正向路径概率和反向路径概率，(\Delta S_{\rm traj}=\ln[P(x)/P_R(x)])。这个量在随机热力学里很核心，模拟上却比平衡态熵麻烦，因为你要记录整条路径的概率，而不只是末态分布。做这类模拟时，随机数质量、路径采样步长、反向路径构造都会影响结果。我的建议是先把平衡态 (g(E)) 和自由能算准，再碰动态量；平衡态都没闭环，动态熵产生很容易变成曲线拟合。

另一个扩展方向是熵产率的粗粒化估计。实际系统自由度太多，不可能全记录，只能选可观测变量做粗粒化，再估计粗粒化熵产生。这时候结果依赖分区方式，不能把它当成唯一物理量。我通常会把粗粒化尺度作为参数扫描，看熵产率是否在某个尺度区间稳定。稳定区间才有解释价值。这个思路和重整化有点像：不是追求一个绝对数，而是看量随分辨率的变化趋势。熵的模拟做到后面，很多时候不是算一个数，而是算一张随参数变化的图。

6.2 我在项目里反复验证的三条经验

第一条经验是：先把“熵是什么”写成一句话，再写代码。热力学熵、冯·诺依曼熵、Rényi 熵、纠缠熵、轨迹熵，每个定义对应的数值对象不同。写代码前把这句话贴在注释里，后面调参时才不会跑偏。我见过最可惜的情况是代码写得很漂亮，采样也收敛，最后发现算的是另一个熵。定义对齐不是形式主义，它决定你该用态密度、密度矩阵还是路径概率。

第二条经验是：所有常数项都要有来源。离散模型用 (2^N) 定标，连续模型用参考温度积分，量子模型用高温极限或基态简并。常数项不确定，绝对熵就不可比。相对熵、比热、自由能差可以绕过常数，但报告里必须写清楚哪些量是绝对的、哪些是相对的。审稿人问起来，你能立刻指出定标步骤，而不是说“整体平移不影响”。

第三条经验是：交叉验证要设计在项目开始，而不是项目结束。先想好小系统基准、独立方法、极限检查，再跑大计算。这样一旦结果不对，你能快速定位是物理效应还是实现问题。熵的模拟尤其如此，因为它不像能量那样有直观的局域图像。多做一次基准，少熬三个晚上。这个习惯听起来朴素，但它是我做过所有数值项目里回报率最高的一条。

最后再分享一个实操小技巧：把每次运行的参数、随机种子、能量窗、最终 (\ln f)、直方图覆盖率写进一个轻量日志文件，和结果图放在一起。过两周回头看，你未必记得当时为什么把窗设成那个范围。日志不需要复杂格式，几行文本就够。熵的模拟结果对参数敏感，留下痕迹，比事后凭记忆重建要可靠得多。