3.4 长程静电与 PME

3.4 长程静电与 PME

本节摘要：库仑势按 1/r 衰减，任何简单截断都会引入可观误差——这是长程静电被单独列为一节的全部理由。本节从"截断为什么是罪"讲起，走一遍 Ewald 分解的思路（不做完整推导，保留可复核的标度论证），讲清 PME 的实空间/倒易空间分工与三个网格参数的配置逻辑，并用一格点间距演算算出 6 nm 盒子的格点规模。读完你应当能解释 rlist、rcoulomb、fourierspacing 三者的关系，并判断一套 PME 配置是否自洽。

3.3 把温压控住了，但静电这块精度短板还悬着。承上：1.2 的六件套里，LJ 项 r⁻¹² 衰减极快、1 nm 外几乎归零，简单截断误差可忽略；静电项 1/r 衰减太慢，1 nm 截断意味着把每个离子对的一半相互作用直接扔掉。启下：本节定下的 coulombtype = PME 会写进 3.5 与 3.6 的所有 mdp，也是 6.3 性能分析里"CPU 倒易负载"的来源。

截断之罪：一个数值实验的教训

把盐水电解质的静电在 0.9 nm 处硬截断会发生什么？历史文献与无数新人的实践给出一致答案：离子对分布函数出现人工峰、介电常数偏离、界面体系的表面张力系统性漂移。罪魁祸首是截断在边界上造成的能量不连续——原子对跨过截断半径的瞬间，相互作用从有限值跳到零，等效于在势能面上凿了一圈悬崖，原子会在悬崖边反复抖动。

补救方案两代：第一代是移位/开关函数（shift/switch），把势能在截断处平滑归零，工程上可用但改变了势能形状本身；第二代就是 PME——不截断，而是把长程部分整体搬到傅里叶空间用网格算，从根上消灭截断误差。现代 GROMACS 只推荐第二代。

Ewald 的拆法：镜像求和的两难如何破

周期性体系中每个电荷面对的是全空间镜像的无穷求和，1/r 项级数收敛极慢。Ewald 的思路是分而治之：给每个点电荷套一团高斯弥散电荷（中和点电荷的奇异性），于是总势能 = "点电荷 + 高斯云"的短程部分 + "高斯云"本身的长程部分 + 自能修正。关键是两部分的收敛速度都可以调：

• 实空间部分：点电荷与高斯云相消，势差按误差补偿函数 erfc(κr)/r 衰减——衰减率由分裂参数 κ 控制，κ 越大实空间收敛越快。GROMACS 里它由 rcoulomb 与相对误差容忍度（verlet-buffer-tolerance 体系）自动定，你不用手填。
• 倒易空间部分：光滑高斯云的和在傅里叶空间里收敛飞快，截断到有限波矢即可。PME（Particle-Mesh Ewald）的工程贡献是：不为每个粒子对求和，而是把电荷"铺"到三维网格上，用快速傅里叶变换（FFT）一次算出整空间的静电势——代价从 O(N²) 降到 O(N log N)。

图 3-3 PME 的空间分工：短程归实空间，长程归网格

演算：格点账与参数自洽检查

fourierspacing = 0.12 nm 的含义是"每个倒易格点负责约 0.12 nm 见方的电荷云"。对 2.3 那个 6 nm 立方盒子：

• 每维格点数 n = 6.0 / 0.12 = 50
• 三维格点总数 = 50³ = 125000
• FFT 代价按 N log N 增长：格点翻倍（间距减半到 0.06）总代价约 2³ × (log 2·50³ 增量) ≈ 8 倍——网格加密的代价比截断加长凶得多

参数自洽检查三问：rcoulomb 是否小于半盒宽（6 nm 盒配 1.0 nm，通过）；pme_order 是否保持 4（六阶只在你真的 chasing 精度且不在乎算力时用）；fourierspacing 是否被 pfme 程序按盒子实际尺寸取整（GROMACS 会把网格调成 FFT 友好的 2、3、5 的倍数组合，log 里打印实际格点数，50 会被取成 48 或 54 附近——看 log 而不是看你的手算）。

还有一条与 3.1 呼应的暗线：rcoulomb 不是孤立参数。Verlet 方案下，实空间误差预算由 verlet-buffer-tolerance 统一掌管，rcoulomb 实际是程序按预算倒推的——你手填 1.0 只是"起点意向"。强行手调（把 tolerance 拉小）时务必同步看 log 里重新报出的 rlist 与格点数，两头的误差要一起缩，只缩一头是白花钱。

何时可以不用 PME

对称的问题也值得回答：什么时候 PME 是浪费？小体系（几百原子、盒子 3 nm 以下）与真空/半真空边界体系，反应场（reaction-field）或纯截断加校正可能更划算——PME 的 FFT 全局同步在小体系上开销占比大。而界面体系（脂膜）用 PME 有个已知的微妙点：周期性在法向同样成立，膜与它的镜像隔着水层"隔空感应"，水层厚度（常用 4 层到 8 层水，约 12 到 22 nm 总盒高）要足以压低这种伪相互作用——7.2 的案例会给出具体的盒高与检验方法。

静电方案的全菜单与取舍

PME 之外的其他 coulombtype 不是不能用，而是各有明码标价：

粗粒化是个重要的例外场景：MARTINI 一族的惯例是反应场加 1.1 到 1.2 nm 截断——粗粒化的静电本来就被大幅平均，反应场的介电修正足够覆盖，而省下的 FFT 开销让上百万原子的介观模拟跑得动。全原子体系则没有这个借口，PME 是默认。

演算：误差预算的分配直觉

PME 的总误差近似为实空间与倒易空间两部分的平方和，两条腿各自能调：实空间误差随 rc 指数下降（由 erfc(κ·rc) 主导），倒易误差随网格密度指数下降（由波矢截断主导）。整体最省算力的配置在两腿误差相等处——一头独大时，另一头的算力全是浪费。GROMACS 用 verlet-buffer-tolerance（默认每原子每步 0.005 kJ/mol）自动配平，手算一遍找感觉：

• 把 tolerance 收紧一个数量级（0.005 → 0.0005），rc 约需从 1.0 增到 1.1 nm（邻居对数涨约三成），fourierspacing 从 0.12 收到 0.10（格点数涨约七成）
• 直觉结论：每要一位精度，两边各付一档代价；只收一边（如只加密网格不涨 rc）省下的误差远比想象少

这也解释了 3.1 的告诫为什么成立：手填 rcoulomb 而不管 tolerance，等于只拧了一边的螺丝。

两个高频疑问

为什么 log 里报的网格数不是 50？ FFT 对 2、3、5 的倍数组合有最快的分解路径，GROMACS 会把 50 这类数取整到邻近的 48 或 54 并在 log 里如实报告——你手算的 50 只是"间距意向"，log 里的数才是事实。核对参数时以 log 为准，这也是"运行前读一遍 log"的习惯之一。

长程色散为什么不用 PME 也行？ LJ 的吸引臂 r⁻⁶ 收敛快，1 nm 外的尾巴用解析尾校正（DispCorr）补即可，误差远小于静电 1/r 的尾巴——这就是 DispCorr = EnerPres 存在的理由。但注意它只修能量与压强的均值，不修涨落，面积压强敏感的脂膜课题要在论文里注明是否开启（7.2 的口径对齐条目）。

速记：LJ 截断无罪、静电截断有罪；PME = 实空间邻居对 + 倒易空间网格 FFT，代价 O(N log N)；三个数——rcoulomb 1.0、spacing 0.12、order 4；格点账用 6/0.12 = 50 起步；精度预算两头配平，log 里的真实格点数才算数。静电这关过了，可以开始走台第一幕：能量最小化。