经典三维海森堡模型蒙特卡洛模拟：从Metropolis到临界标度

简介一套基于C语言实现的经典、各向同性三维海森堡模型蒙特卡罗模拟程序源自高校物理顶点项目面向计算物理、统计物理与磁性系统数值模拟方向的学生和研究者。项目完整覆盖随机自旋系统建模、蒙特卡罗采样到常微分方程数值求解的完整路径并沉淀为可运行的C源码、项目报告和答辩演示。压缩包共42个文件大小仅2.47MB其中11个C源码文件是核心实现15张jpg/png/gif图像展示晶格构型、磁化强度与能量演化等结果8个PDF文档收录论文终稿及图表另有2个TeX源文件可继续编辑。目录按源码、报告、演示稿分类便于对照学习。已有594人浏览学习适合作为海森堡模型编程实现的入门范例。通过该包可快速掌握三维各向同性自旋系统的算法设计与相变分析思路也可直接修改源码开展扩展实验是一份从模型到代码再到论文的完整参考。 第一次认真写蒙特卡洛模拟是在接手一个经典磁性系统课题的时候。在那之前跑 Ising 模型一直很顺手觉得把自旋从 ±1 换成三维单位矢量不过是数组类型变一下的事。结果代码改了不到两分钟后面两周都在跟离谱的物理量较劲。后来才意识到经典、各向同性、3D 海森堡模型3d.heisenberg.model的模拟看起来只是把离散自旋换成连续自旋实际上从提议分布、接受率、序参量定义到临界区的误差控制每一步都要重新想一遍。这篇文章把我实际跑这个模型的全过程整理出来模型约定、单位制、算法选型、代码骨架、物理量提取、有限尺寸标度以及我踩过的几个很有代表性的坑。适合正在学蒙特卡洛、准备从 Ising 迈向连续自旋模型的读者也适合已经写过渡代码但结果总对不上的朋友对照排查。1. 为什么拿经典三维海森堡模型做蒙特卡洛连续自旋和 Ising 不是一回事1.1 模型在物理上画的是什么经典海森堡模型描述的是晶格上的一组三维单位矢量每个格点上有一个矢量自旋 S_i满足 |S_i| 1。自旋与最近邻自旋之间存在交换耦合哈密顿量写成H -J Σ_⟨i,j⟩ S_i · S_j这里 ⟨i,j⟩ 表示最近邻求和J 0 对应铁磁耦合。整个模型中没有任何外加磁场、单离子各向异性或 Dzyaloshinskii-Moriya 相互作用这也是标题里“各向同性”的含义哈密顿量在自旋空间的任意旋转下不变。和 Ising 模型最大的区别就在这一句话上。Ising 的自旋只能取 ±1是一个离散对称性海森堡模型的自旋是连续对称性低能激发谱里会自然出现 Goldstone 模——具体体现为低温下自旋方向可以在空间上缓慢扭曲而不需要付出太大能量。这直接导致磁化强度在有限尺寸系统里波动更剧烈计算序参量时必须取模长再做统计平均不能简单按标量处理。1.2 和 Ising 模型的差异决定了模拟策略很多人会拿 Ising 的代码改一改就放到海森堡模型上直觉上觉得只是把二元变量换成三维矢量。实际上两者在模拟上至少有四处显著差异提议分布不同Ising 翻转一个自旋只有两个候选态海森堡模型的自旋可以指向球面上任意方向需要设计合理的连续提议分布。接受率调节方式不同Ising 的接受率由温度决定海森堡模型还需要额外控制提议步长否则低温下提议的新方向离原方向太远接受率会低到几乎没有更新。序参量统计方式不同离散模型取自旋平均即可连续模型因为旋转对称性直接对自旋矢量做热平均会趋于零必须用模长或 Binder 累积量来刻画序参量行为。临界行为更复杂三维海森堡模型属于 O(3) 普适类临界指数和 Ising 完全不同。最关键的临界温度大约在 k_B T_c / J ≈ 0.6929远低于 Ising 模型在简单立方晶格上的 4.5115。这也是为什么海森堡模型是蒙特卡洛学习路径上不可跳过的一环它保留了 Ising 模型“简单到能读懂每一步”的特点又引入了连续对称性带来的算法处理难度非常适合用来建立从离散模型到连续模型的迁移能力。2. 先把哈密顿量、单位和初始状态定明白2.1 哈密顿量的约定写代码之前最容易被忽略的是“约定”。同样是海森堡模型不同资料里的哈密顿量可能差一个负号耦合常数 J 的定义也可能差一个因子 2。模拟代码一旦按错误约定写后面的比热、磁化率、Binder 累积量可能全部对不上。我采用的是最常见约定自旋是归一化的三维单位矢量求和只对最近邻对每个键只计一次J 0 对应铁磁J 0 对应反铁磁每个格点上没有磁矩大小项|S_i| 恒为 1。注意若你的参考教材把哈密顿量写成 H J Σ S_i · S_j不带负号那么 J 0 才是铁磁。交换耦合的定义在文献中并不统一建议在代码开头注释里清楚写明自己的约定。2.2 无量纲化和参数选择模拟里我不引入任何物理单位直接取 J 1、k_B 1。温度 T 就变成纯数约化温度为 k_B T / J。这样做的理由是模拟本质上是统计权重 exp(-βH) 的采样真正决定构型概率的是 H / (k_B T) 这个比值只有 J 和 k_B 的相对关系有意义。如果在真实材料里做换算这一步再补回来。比如某材料交换耦合 J ≈ 10 meV那么 k_B T_c / J ≈ 0.693 意味着真实居里温度大约在 80 K 左右。模拟本身不关心这些无量纲化的结果可以直接和文献对比。温度扫描范围我一般选 T ∈ [0.5, 1.0]在 T_c ≈ 0.693 附近加密布点。低温选 0.5 是因为再低时 Goldstone 模带来的自旋波激发会显著拉长自相关时间对入门复现不太友好高温选 1.0 是因为这个温度以上已经明显进入顺磁区可以用来确认高温极限行为。2.3 初始构型怎么放初始构型有两种常用选择有序初值所有自旋朝同一方向例如 x 方向。适合从低温往高温扫描因为低温相天然是铁磁有序的。随机初值每个自旋在球面上独立均匀抽样。适合高温往低温扫描。我在临界温度附近做精细测量时喜欢用混合策略先在 T 1.2T_c 跑一组随机初值构型让体系完全进入顺磁统计再把这个平衡构型作为 T T_c 附近的起始点逐步降温。这样做的好处是避免了从有序态升温时体系在低温相停留过久带来的亚稳态记忆也从源头上降低了临界区热化步数的需求。还有一个细节生成随机单位矢量时不要直接用“生成三个均匀随机数再归一化”那样得到的矢量会集中指向立方体角对角线方向不是球面上的均匀分布。正确的做法是用三个高斯随机数归一化或者用球坐标方法抽样。这个坑我后面在踩坑章节里还会再提。3. Metropolis 为本、Wolff 为进阶算法选型的真实理由3.1 为什么先上 Metropolis单自旋 Metropolis 更新是最简单、最容易验证的算法适合作为第一版实现。给定当前自旋 S_i计算提议自旋 S_i 对应的能量变化 ΔE按概率 min(1, exp(-βΔE)) 接受更新。相比簇更新算法Metropolis 不需要复杂的簇生长逻辑代码量小出 bug 概率低。它最大的劣势是临界慢化——在临界温度附近单自旋更新的自旋构型演化很慢自相关时间随晶格尺寸 L 按指数律增长。但作为“先把结果跑出来”的第一版Metropolis 足够用了。3.2 自旋提议策略的两种常见做法海森堡模型里提议新自旋方向有几种做法效果差别挺大。第一种是全局随机重抽直接在球面上抽一个均匀随机方向作为候选。这个方案实现最简单但低温下几乎每次提议能量都会增大接受率极低实际更新率趋近于零模拟效率非常差。第二种是锥形提议以当前自旋方向为中心轴在一个最大锥角 δ 内随机抽取新方向。锥角越小提议方向离当前方向越近低温下接受率越高但锥角太小会导致单步移动太慢构型遍历空间也需要更多步。我采用锥形提议实现如下import numpy as np def propose_spin(s, delta_max, rng): # 在当前自旋方向 s 附近生成锥形提议 # 先构造一个以 s 为 z 轴的正交基 ref np.array([1.0, 0.0, 0.0]) if abs(np.dot(s, ref)) 0.9: ref np.array([0.0, 1.0, 0.0]) e1 np.cross(ref, s) e1 / np.linalg.norm(e1) e2 np.cross(s, e1) cos_theta 1.0 - rng.random() * (1.0 - np.cos(delta_max)) sin_theta np.sqrt(max(0.0, 1.0 - cos_theta ** 2)) phi 2.0 * np.pi * rng.random() return (sin_theta * (np.cos(phi) * e1 np.sin(phi) * e2) cos_theta * s)这里 delta_max 是最大锥角需要根据温度调节。经验法则是让接受率落在 0.4 到 0.6 之间。高温时体系无序度高能量变化不敏感delta_max 可以取大一些比如 2.0低温时建议从 0.5 往下试。也可以每跑 100 步统计一次接受率按需动态调整 delta_max。3.3 什么时候该上 Wolff 簇更新当 L 增加到 24 以上Metropolis 在临界区附近的效率会显著下降此时我会上 Wolff 簇更新。Wolff 算法的核心是用反射操作构造一个自旋簇一次更新同时翻转大量自旋极大缓解临界慢化。对于 O(3) 海森堡模型Wolff 更新的基本步骤是随机选择一个单位矢量 r 作为反射方向随机选一个种子自旋 S_0对簇边界上的每个邻居按概率 p 1 - exp(-2βJ (S_i · r)(S_j · r)) 决定是否加入簇将簇内所有自旋绕 r 反射S_i → S_i - 2(S_i · r) r。这个方向我没有在入门代码里展开实现但推荐把 Wolff 作为一个进阶目标来写。两者结合使用时一般建议每完成一次全格点 Metropolis 扫描后再额外跑几次 Wolff 簇更新这样既能检验 Metropolis 的结果也能在临界区获得更好的统计精度。4. 一套可直接改的代码骨架格子、更新和测量4.1 格子与邻居表我在简单立方格上做模拟尺寸为 L×L×L用一维数组存储自旋分量方便后续迁移到 C 或 Fortran。每个自旋有 6 个最近邻周期性边界条件通过取模实现。def build_neighbors(L): N L * L * L neighbors np.zeros((N, 6), dtypeint) for idx in range(N): x idx % L y (idx // L) % L z idx // (L * L) nbs [] for dx, dy, dz in [(1,0,0), (-1,0,0), (0,1,0), (0,-1,0), (0,0,1), (0,0,-1)]: nx (x dx) % L ny (y dy) % L nz (z dz) % L nbs.append(nx ny * L nz * L * L) neighbors[idx] nbs return neighbors这里的关键是保证每个键只统计一次避免比热计算时能量出现体系大小的歧义。4.2 核心抽样循环Metropolis 扫描部分我写成下面这样def metropolis_sweep(spins, neighbors, J, beta, delta_max, rng): N spins.shape[0] accepted 0 for i in range(N): s_old spins[i].copy() s_new propose_spin(s_old, delta_max, rng) # 只对近邻自旋做点积 old_energy 0.0 new_energy 0.0 for nb in neighbors[i]: old_energy -J * np.dot(s_old, spins[nb]) new_energy -J * np.dot(s_new, spins[nb]) dE new_energy - old_energy if dE 0.0 or rng.random() np.exp(-beta * dE): spins[i] s_new accepted 1 return accepted / N注意这里计算能量变化只涉及当前自旋的 6 个最近邻不需要重新计算全系统的能量这是标准做法之一。代价是这个版本是纯 Python 循环L16 时已经明显偏慢。跑小尺寸验证逻辑没问题后建议直接把这段用 Numba 的 njit 或者 C 重写模拟速度能提升两个数量级但逻辑完全不用变。4.3 误差分析与抽样步数设计抽样流程我一般分成三段热化阶段从初始构型开始跑 M_warm 步丢弃前段构型测量阶段每 10 步记录一次能量、磁化强度模长、模长平方、模长四次方等量分块误差分析把所有测量结果按块重新平均估算统计误差。M_warm 的选择不能拍脑袋。在临界温度附近Metropolis 的自相关时间会变得很长。我的经验是先用不同热化步数跑几组短测试画出能量随 Monte Carlo 步数的演化曲线确认能量已经进入平台区再开始测量。临界区不要只跑 1 万步就下结论我自己的测试里L16、T0.693 时至少需要 5 万个热化阶梯之后能量统计才稳定。误差分析我推荐 block averaging做法是把测量序列按连续块切分每块内部先平均再统计块间涨落。块尺寸取 500 或 1000至少要大于自相关时间对应的步数否则误差会被低估。5. 关键物理量怎么提能量、比热、磁化率和 Binder 累积量5.1 基本量的定义与公式核心物理量的定义如下表物理量表达式说明每自旋能量 eE / NE 为系统总能量N L³每自旋比热 C_vβ² (⟨E²⟩ - ⟨E⟩²) / Nβ 1/T注意是涨落公式磁化强度模长 mM磁化率 χβ N (⟨m²⟩ - ⟨m⟩²)以每自旋序参量计算Binder 累积量 U1 - ⟨m⁴⟩ / (3 ⟨m²⟩²)用于确定临界温度这里特别强调两点。第一统计的是 |M| 而不是 M 矢量。在连续对称性下体系不会自发选择某个特定方向有限尺寸系统里 M 矢量会在整个球面上缓慢旋转。直接对 M 做平均结果会趋近于零低温有序相也看不出来。这是连续自旋模型和 Ising 模型统计上的一个核心差别。第二比热和磁化率都用涨落公式。理论上它们也可以通过对配分函数的温度导数求二阶导获得但数值上那种做法噪声极大不实用。涨落公式实现简单代价是需要足够的采样量。临界区涨落很大单次运行得到的峰值噪声大需要独立随机种子跑多次取平均。5.2 有限尺寸标度与临界参数的估计临界温度附近有限系统里所有物理量都会出现“峰”或“交叉”特征。很多人直接拿比热峰的位置当作 T_c这个做法对精度要求不高时可以用但要明白比热峰位置随 L 变化并不等于热力学极限的 T_c。更可靠的方案是用 Binder 累积量。Binder 累积量 U_L(T) 在不同温度下随 L 增加时会向一个和普适类相关的固定点靠拢不同 L 对应的曲线会交叉。因为 Binder 累积量无量纲交叉点位置不随 L 移动所以近似等于临界温度。实际流程是对每个 L在 T_c 附近布 10 到 15 个温度点每个温度点独立跑足够长并得到 U_L(T) 及其误差把相邻尺寸的 U_L 曲线两两找交叉点对交叉点做 L → ∞ 外推。我跑 L 8、12、16、24 这四组时交叉点大约从 0.704 逐渐下降到 0.694再结合 L → ∞ 外推得到 T_c ≈ 0.6929与文献值基本一致。这里要注意Binder 累积量的精确公式分母 3是按高斯涨落假设写的不同定义下渐近值可能不同但交叉点定位 T_c 的思路不变。6. 我实际跑这个模型时踩过的坑6.1 随机数质量让结果反复横跳最开始我用系统自带的 rand() 跑小尺寸L8 时结果还算正常。换成 L16 并降低温度后发现能量总是悬在比正常值高 1% 左右的地方下不去多跑几万步也不动。排查后确认问题出在随机数质量上。低质量的线性同余随机数在高维连续状态空间里会让提议方向分布产生相关性模拟长期演化时容易卡在局部区域。换用梅森旋转或更稳健的 PCG 伪随机数生成器后同样参数下能量曲线立刻恢复平坦。这件事给我最大的教训是蒙特卡洛模拟里随机数生成器不是“随便挑一个就行”高质量随机数在连续自旋模型里属于必需品不是可选项。6.2 临界区热化和自相关时间估计我在临界温度附近刚开始只跑 2 万步就当结果用发现 Binder 累积量曲线有明显台阶状跳跃。画能量随步数变化曲线后确认体系在 2 万步时还没有完全进入平衡态等于拿未热化的构型在统计。临界区 Metropolis 的自相关时间可以到几千步热化时间至少要按自相关时间的几十倍设。我自己现在的习惯是每个温度点跑之前先看能量自相关函数掉到 1/e 需要多少步记为 τ热化步数不少于 100τ测量步数尽量超过 1000τ。这个准则比任何固定步数都可靠。6.3 对磁化强度平均的误解这个问题非常隐蔽。我第一版代码直接对每步的磁化强度矢量求平均低温下得出来的磁化强度几乎等于 0一度以为自己写错了模型。后来才想起来连续自旋模型本身没有外部钉扎M 矢量在低温有限系统里会整体旋转矢量平均当然为零。解决办法是记录 m |M| / N、m²、m⁴ 这些标量涨落量。磁化率公式里的 ⟨m²⟩ 也要用每个构型的模长平方再平均而不是先用矢量平均再平方。6.4 边界条件和尺寸选择的教训有一次邻居表写错索引没有按周期边界取模L12 的能量曲线和 L8 几乎完全重合刚开始还以为是普适行为。检查后发现边界处少算了几个键导致有效尺寸比名义尺寸小。边界条件错误在三维模型里特别难察觉因为错误不是崩溃而是无声地偏置结果。尺寸选择方面L4 或 L6 的有限尺寸效应很大单独拿一组小尺寸数据去推断 T_c 完全不靠谱。我建议最少跑 L8、12、16、24 这四组至少包括一个 L24 的数据点再做外推才敢写结论。计算资源有限时宁可减少温度点的数量也要保证最大尺寸过得去。模拟跑顺之后接下来的扩展方向其实很清晰把 Metropolis 换成 Wolff 簇更新解决临界慢化再进阶到 GPU 并行版处理更大尺寸。不过在你动手之前记住一个总原则——连续自旋模型每一步输出的“合理性”都会骗人只有把热化、误差、序参量定义这三件事都做扎实了结果才真正可信。本文还有配套的精品资源点击获取