3.3 温度与压力耦合对照
文档摘要
3.3 温度与压力耦合对照 本节摘要：恒温器与恒压器决定你的模拟名义上处于哪个统计系综：不是所有算法都产生正确的正则或等压系综，有些只适合"往目标态赶路"，有些才能"驻留"。本节用一张对照表把三类恒温器与三类恒压器的机理、可逆性与涨落保真度讲透，用涨落公式演算热容，给出 taut、taup 与压缩率的推荐值及其来历，并立下"平衡期可以松、生产期必须严"的阶段纪律。 3.2 定下了步长与积分器，但一个只有牛顿方程的模拟是孤立系统（NVE 系综）：总能量严格守恒，温度只会围绕初值漂移。真实实验在恒温槽与大气压里做，所以模拟也要把体系"泡进"虚拟恒温槽与虚拟活塞——这就是耦合算法的全部职责。选择算法的核心判据只有两条：系综对不对（统计力学严格性）与涨落真不真（物理可观测量是否失真）。
3.3 温度与压力耦合对照
本节摘要：恒温器与恒压器决定你的模拟名义上处于哪个统计系综：不是所有算法都产生正确的正则或等压系综，有些只适合"往目标态赶路"，有些才能"驻留"。本节用一张对照表把三类恒温器与三类恒压器的机理、可逆性与涨落保真度讲透，用涨落公式演算热容，给出 tau_t、tau_p 与压缩率的推荐值及其来历，并立下"平衡期可以松、生产期必须严"的阶段纪律。
3.2 定下了步长与积分器，但一个只有牛顿方程的模拟是孤立系统（NVE 系综）：总能量严格守恒，温度只会围绕初值漂移。真实实验在恒温槽与大气压里做，所以模拟也要把体系"泡进"虚拟恒温槽与虚拟活塞——这就是耦合算法的全部职责。选择算法的核心判据只有两条：系综对不对（统计力学严格性）与涨落真不真（物理可观测量是否失真）。
恒温器三家：松的、严的、介于其间的
| 算法 | 机理一句话 | 系综严格性 | 涨落保真 | 建议用途 |
|---|---|---|---|---|
| Berendsen | 速度按偏差比例重标定，指数趋近目标 | 不严格 | 压低（约 1/2 量级） | 仅限平衡初段"灭火" |
| V-rescale | 重标定加随机项，来自正则系综的严格随机动力学 | 严格（正则） | 基本正确 | 平衡与生产通吃，GROMACS 默认推荐 |
| Nosé-Hoover | 引入一个"热浴动量"变量的确定性反馈 | 严格（正则） | 正确 | 生产可用；需配 tau 与链长调好 |
| (参考) 不耦合 | 无 | NVE | — | 检验积分器守恒性时用 |
Berendsen 的问题值得多讲一句：它把动能涨落强行压平，而动能涨落本身就是物理——C_V 就藏在涨落里（见下面演算）。用 Berendsen 跑生产，你测的热力学量系统性偏差。历史上大量老论文用它，复现时换成 V-rescale 往往连观测量都变，这不是 bug 是修正。
tau_t 怎么取：默认推荐 0.5 ps。取值的权衡是——太小（如 0.01 ps）接近每步硬扳，变成另一种 Berendsen，涨落被压、局部过阻尼；太大（如 5 ps）控温形同虚设，体系长时间偏离目标温度。经验法则：tau_t 取体系本征弛豫时间的几分之一即可，溶剂主导的体系 0.5 ps 几乎永远合适。
恒压器三家：盒子的呼吸方式
| 算法 | 机理一句话 | 系综严格性 | 配套建议 |
|---|---|---|---|
| Berendsen | 盒子按压差比例缩放，稳但假 | 不严格 | 同名恒温器的问题同源，只配平衡前段 |
| C-rescale | 胞胞随机压控，严格等压系综 | 严格 | 平衡与生产通吃，新默认 |
| Parrinello-Rahman | 盒矢量的动力学方程，含惯性振荡 | 严格（NPT） | 生产首选之一；平衡初期易震 |
| (各向异性) ref-p 变量 | 每方向独立压强或压缩率 | 严格 | 脂膜（面内 vs 法向）专用，7.2 详述 |
恒压的 tau_p 推荐值 5 ps（C-rescale 可到 2 ps），比 tau_t 大一个量级，因为盒形变量是集体慢变量，控太快会把声波硬生生按进盒子。压缩率 compressibility 填 4.5e-5 bar⁻¹——这就是液态水的等温压缩率；对聚合物熔体可取 1e-4 到 5e-5 之间，它只影响耦合器"信不信压差"的增益，不改变目标态本身。
图 3-2 温压耦合对照矩阵：一张图记住六种选择
演算：热容藏在能量涨落里
正则系综里，势能涨落与定容热容之间有严格关系（涨落耗散的又一实例）：
C_V = (⟨U²⟩ − ⟨U⟩²) / (k_B T²)
用 03-1 那套 PEO 溶液体系的假想 edr 数据演算：350 K 下势能均值 ⟨U⟩ = −1.286e5 kJ/mol，标准差 σ_U = 410 kJ/mol（gmx energy 输出的 RMSD 栏就是它），则
C_V = 410² / (0.008314 × 350²) ≈ 168100 / 1018.5 ≈ 165 kJ/(mol·K)
对约 3500 个有效自由度的体系，经典均分预言 C_V = N_df · k_B / 2 = 3500 × 0.008314 / 2 ≈ 14.6 kJ/(mol·K)——两值差一个量级？先别慌：edr 的势能涨落里混进了约束反弹与耦合算法的贡献，且单位是每摩尔模拟体系而非每摩尔物质。这笔账真正想教的是方法论：任何"从涨落算热力学量"的流程，第一步都是核对系综正确性——如果 3.1 用了 Berendsen，这 165 就是个系统性缩小的假数。涨落保真不是洁癖，是定量结论的地基。
阶段纪律与易错点
把两节内容合成一条作业流水线：平衡初段（前 50 ps 量级）可以用 Berendsen 温压组合快速灭火；中段换 V-rescale 继续 NVT/NPT；正式生产用 V-rescale 加 C-rescale 或 Parrinello-Rahman。换耦合算法不需要重新平衡太久，但要留出新的收敛观察窗（3.6 的判据照常跑）。
三个易错点：其一，tau_t 填了一个值却分了两个组（组数与 tau/ref 列数必须相等，grompp 会报错拦截）；其二，脂膜体系用各向同性恒压——双层法向被压塌、面积压强失真，必须 ref-p = anisotropic 并把面内压缩率设小；其三，从别人的旧 mdp 继承了 Berendsen 跑生产，数据看起来"挺稳"其实涨落全是假的——这种错误不会崩，只会让结论慢性中毒，比崩溃更可怕。
tau 参数的量级表：抄作业要抄对场景
| 体系 | tau_t 推荐 | tau_p 推荐 | 说明 |
|---|---|---|---|
| 全原子溶剂体系 | 0.5 ps | 5 ps | 行业默认，几乎通用 |
| 联合原子/粗粒化熔体 | 0.5 到 1.0 ps | 5 到 10 ps | 集体弛豫更慢，恒压放宽 |
| 脂膜双层 | 0.5 到 1.0 ps | 面内 5 ps、法向 15 ps | 各向异性分开（7.2） |
| 气液界面/真空占比大 | 0.5 ps | 10 ps 以上 | 低密度区反馈噪声大 |
| NVE 检验运行 | 不耦合 | 不耦合 | 只为验证积分器守恒性 |
这张表的用法是"起点"而非"结论"：tau 的最终取值要看 3.6 判稳脚本的输出——温度均值贴不上目标就调小 tau_t，压强反馈让盒子抖得厉害就放宽 tau_p。参数与判据联动，才是调参而不是抄参。
两个高频疑问
为什么平衡期也建议 V-rescale 而不是 Berendsen 一路到底？ 因为平衡期同样在统计——你判断"平衡了没"用的温度、密度均值，全靠涨落是真实的。Berendsen 压平涨落后，3.6 那套判稳脚本会给出过于乐观的"平稳"假象，可能把没平衡的体系放行。V-rescale 两头兼顾：收敛速度接近 Berendsen（随机项帮助翻越动能垒），统计性质又严格，所以它是"从平衡一路用到生产"的省心选项。
表面张力耦合怎么写？适合谁？ pcoupltype = surface-tension，给面内一个目标张力（通常 0，让膜自己找面积）、法向控压强。它适合"想让膜自由决定面积"的课题；但要注意报告面积压强时必须注明耦合方式——同一张膜在 surface-tension 与 anisotropic 下平衡面积可以差出几个百分点，这不是谁对谁错，是系综不同。
速记：恒温器选 V-rescale 起步，热浴动量一族备选，Berendsen 只配灭火；恒压器选 C-rescale 起步，P-R 做生产，各向异性专治脂膜；tau_t 0.5 ps、tau_p 5 ps、水压缩率 4.5e-5 是三个可以直接抄的数；涨落是物理，压平涨落等于删数据。控住了温压，下一节解决静电这个精度大坑。