1.2 势能函数与力场的数学骨架

1.2 势能函数与力场的数学骨架

本节摘要：力场是一套"函数形式 + 参数表"的组合：函数形式规定能量怎么随几何变化，参数表给出每种原子组合的具体数值。骨架由六大件构成——键长、键角、二面角、非正常二面角四项成键项，加 Lennard-Jones 与静电两项非键项。本节逐项拆解每一件在守什么、参数量纲是什么，并用一段 Python 手算一个乙烷片段的总能量，让你在进入第2章选型之前先摸过参数本身。

上一节把力定义成势能的负梯度，但那张势能表内部长什么样还没有拆开。这非常关键：你以后在 mdp 文件里调的每一个截断、在拓扑里看到的每一条 bond 定义，都是这副骨架上的关节。本节从分子的几何讲起，按"成键四件套 + 非键两兄弟"的顺序拆，最后把六项加总做一次真实演算。拆完这节，2.1 的力场选型就不再是背名词，而是比较不同厂家给同一副骨架配的不同参数。

全景：分子是一根挂着电荷的弹簧床

把分子摊开看，几何自由度只有三类：两原子间的距离（键长）、三原子构成的夹角（键角）、四原子构成的扭转（二面角）。力场给每一类几何配一项能量函数，再对分子内所有非键原子对配上排斥—吸引与静电两项。一个典型力场写出来是：

U_total = Σ_bonds k_b(r − r₀)² + Σ_angles k_θ(θ − θ₀)² + Σ_dihedrals k_φ[1 + cos(nφ − φ₀)] + Σ_impropers k_ψ(ψ − ψ₀)² + Σ_ij [ 4ε_ij( (σ_ij/r)¹² − (σ_ij/r)⁶ ) + q_i q_j / (4πε₀ r) ]

前四项只对拓扑里连了线的原子生效，后两项对"隔了三个键以上"的原子对生效（1-2、1-3 通常被排除，1-4 对按力场特殊处理）。这个"谁能对谁作用"的规则表叫排除表，是拓扑文件里看不见却致命的部分——2.2 生成拓扑时它会自动写好，但你得知道它在哪。

图 1-2 势能六大项与它们守住的几何

成键四件套：几何的三层守卫

键长项与键角项都是简谐形式，物理含义最直白：偏离平衡位置越远、回弹力越大。参数量纲要盯紧——键长力常数是 kJ/(mol·nm²)，角度力常数是 kJ/(mol·rad²)。演算一笔：聚环氧乙烷（PEO）里醚氧两侧的 C–O–C 键角，典型参数 θ₀ ≈ 1.87 rad（约 107°）、k_θ ≈ 400 kJ/(mol·rad²)。若某一步积分把角度压到 95°（1.658 rad），这一项的能量是

U = 400 × (1.658 − 1.87)² = 400 × 0.0449 ≈ 18.0 kJ/mol

对照室温热能 k_B T ≈ 2.48 kJ/mol，这一个角的形变就吃掉约 7 倍热能——能量最小化没做好的体系里，这种"高能关节"成片存在，正是 4.3 里"原子飞出盒子"类崩溃的第一来源。

二面角项是四件套里唯一的多谷形式：cos 展开允许同时叠加多个周期（比如 n=3 主峰加 n=1、n=2 修正），这正是它能同时刻画"反式更稳、旁氏略高、左右旁氏对称"这种构象偏好的原因。PEO 链的柔性、聚乙烯链的锯齿构象，全由这一项的参数组合决定。非正常二面角则相反，它通常只有一个平衡值，用简谐形式把芳香环死死按在平面上，或把手性中心的构型锁住——参数失配时环会"翻转"，你在 7.2 的脂膜案例里会见到胆固醇环被压塌的实例。

非键两兄弟：排斥吸引与静电

Lennard-Jones 项的 12 次幂臂负责泡利排斥，6 次幂臂负责色散吸引；σ 是"两原子刚好互不理睬"的距离，ε 是阱深。两种原子的交叉参数通常不直接给，而是用混合规则算出。GROMACS 默认的几何混合：

σ_ij = √(σ_i · σ_j)，ε_ij = √(ε_i · ε_j)

演算：PEO 醚氧（σ=0.3031 nm，ε=0.6776 kJ/mol，取 OPLS-UA 风格数值）与水氧（σ=0.3166 nm，ε=0.6502 kJ/mol）交叉：

• σ_OC = √(0.3031 × 0.3166) ≈ 0.3098 nm
• ε_OC = √(0.6776 × 0.6502) ≈ 0.6641 kJ/mol

注意另一种常见规则 Lorentz-Berthelot 用算术平均混 σ：σ_ij = (σ_i+σ_j)/2 = 0.3099 nm——差别很小，但 ε 在某些力场里也用算术平均，结果可差百分之几十。混用什么规则必须与力场配套，混错了等于换了一个力场，这是 2.1 参数核对清单里的固定栏目。

静电项最朴素也最麻烦：1/r 衰减太慢，任何截断都会引入可观的误差，这一条直接催生了 3.4 的 PME 方法。这里只需记住：力场给每个原子分配一个固定部分电荷（fixed charge），它不随构象变化——极化效应被平均摊进了参数。对水、脂膜这类强极化体系，这正是力场系统误差的主要来源之一。

函数形式不止一套：三大门派的骨架差异

同样的六件套，写法有门派之分，这也是 2.1 选型表的底层逻辑：

门派差异看似只在符号层面，实际会改变构象平衡：同一条十二烷在 AMBER 与 CHARMM 下的反式比例可差好几个百分点。结论只有一条——整个体系从头到尾用同一门派，不要拿 CHARMM 的脂膜配 AMBER 的水，除非你明确知道自己在做什么（7.2 的脂膜案例会演示正确的混搭姿势：全 CHARMM 化）。

演算：Python 手算一个乙烷片段的能量

把六项全部加总是什么体验？下面脚本对一个简化乙烷（3 个重原子模型：C–C 键 + 一个 H–C–C 角 + 一个 C–C–H–H 二面角 + 两碳间 LJ + 两电荷静电）算总势能。所有参数为演示值，量纲与 GROMACS 一致：

import numpy as np

KB = 0.0083144621  # kJ/(mol·K)

# ---- 成键参数（演示量级，真实力场查对应 ff 目录的 *.itp）----
kb, r0   = 280000.0, 0.153     # C-C 键：kJ/mol/nm2, nm
kth, th0 = 305.0, 1.911        # H-C-C 角：kJ/mol/rad2, rad
kphi, n, phi0 = 0.650, 3, 0.0  # 二面角：kJ/mol, 周期数, 相位 rad

# ---- 非键参数 ----
sigma_cc, eps_cc = 0.340, 0.36   # C 与 C 的 LJ
q1, q2 = 0.05, -0.05             # 部分电荷 e
EPS0 = 0.138935456               # kJ·nm/(mol·e^2)，GROMACS 用的库仑常数

def u_total(r, theta, phi, d_cc=0.153):
    u_bond = kb * (d_cc - r0) ** 2
    u_ang  = kth * (theta - th0) ** 2
    u_dih  = kphi * (1 + np.cos(n * phi - phi0))
    u_lj   = 4 * eps_cc * ((sigma_cc / d_cc) ** 12 - (sigma_cc / d_cc) ** 6)
    u_coul = EPS0 * q1 * q2 / d_cc
    return u_bond, u_ang, u_dih, u_lj, u_coul

# 平衡构象附近小幅形变，逐项打印
for d_cc, theta, phi in [(0.153, 1.911, 0.0),   # 平衡
                          (0.163, 1.911, 0.0),   # 键拉伸 0.01 nm
                          (0.153, 1.750, 0.0),   # 角压缩约 9 度
                          (0.153, 1.911, np.pi/3)]:  # 二面角扭 60 度
    terms = u_total(None, theta, phi, d_cc)
    print(f"d={d_cc:.3f} th={theta:.3f} phi={phi:.2f} -> "
          f"bond {terms[0]:8.3f} ang {terms[1]:7.3f} dih {terms[2]:5.3f} "
          f"lj {terms[3]:7.4f} coul {terms[4]:7.4f} | 总 {sum(terms):9.3f} kJ/mol")
​

对照输出可以看清三件事：键拉伸 0.01 nm 就贡献约 28 kJ/mol（键是最硬的关节）；角压 9 度约 10.4 kJ/mol；二面角转 60 度不足 1 kJ/mol——软关节先动、硬关节后动，这正是能量最小化（3.5）总是先松二面角、最后才动键长的物理原因，也是步长受键振动（而非二面角）限制的原因。

易错点与本节速记

三个易错点：量纲张冠李戴（键长力常数误当 kJ/mol 用，能量会差三个量级）；混用不同门派的混合规则与 1-4 处理；忽略排除表导致双计数。前两个本节已讲，第三个在 2.2 的拓扑实操里用 pdb2gmx 的输出当面验证。

速记：六件套 = 四成键（键长、角、二面角、非正常二面角）+ 两非键（LJ、静电）；键硬角中二面角软；交叉 LJ 参数用几何规则混，但必须与力场配套；部分电荷不随构象变化是经典力场的核心近似。带着"参数落在哪个文件、谁负责核对"这个问题，进入 2.1。