全原子分子动力学模拟原理：力场选择、时间步长与系综耦合的物理账本

全原子分子动力学模拟原理——把牛顿运动方程用飞秒级步长做数值积分，通过统计力学连接微观轨迹和宏观可观测量——这套理论框架本身是简洁的。但在实际模拟中，三个基础决策环环相扣：力场定义势能面、时间步长定义相空间采样精度、系综耦合定义温度和压力的物理意义。任何一个决策出问题，轨迹上的累积偏差会越来越大，最终让统计结果偏离物理真实。

力场选择：可极化还是固定电荷

全原子力场按极化处理分两类：固定电荷力场（AMBER ff14SB、CHARMM36、OPLS-AA）和可极化力场（AMOEBA、Drude振子）。固定电荷力场把每个原子的部分电荷固定在一个数值上，计算量小，但对静电环境变化敏感的体系（如蛋白内部强极化的活性位点，或离子通道中的钠钾选择性）会漏掉局域极化效应。

团队做过一个水合Ca²⁺离子的案例。CHARMM36力场下，Ca²⁺第一水合层的径向分布函数在2.4 Å和4.6 Å出现两个主峰，配位数~7.2，与EXAFS实验值6-8吻合。但如果用TIP3P水模型配合固定电荷的Ca²⁺，第一峰的精确位置和峰高对Lennard-Jones参数σ和ε非常敏感。这也是为什么力场开发者要反复迭代非键参数——固定电荷的局限性是通过参数拟合来弥补的，而不是物理上被消除。

可极化力场的优势在离子选择性领域最明显。KcsA钾通道中K⁺/Na⁺选择性比高达10000:1——固定电荷力场可以大致重现选择性，但具体的能量差异（~5 kcal/mol）只有在可极化框架下才算得准。团队在是否升级到Drude力场这件事上踌躇过——计算量增加5-8倍，但离子通道项目的结论可靠性对药理推断有直接影响。最终选择了可极化力场，用GPU集群跑了200 ns，得到的选择性自由能profile与实验吻合度明显优于固定电荷结果。

时间步长：2 fs是黄金标准吗

全原子分子动力学模拟原理的数值积分多数用Leap-Frog或Velocity Verlet算法，时间步长设2 fs是标准推荐——因为氢原子X-H键的振动周期约10 fs，2 fs的步长可以稳定积分。去除非键截断外的长程静电（PME）每步都算一次，对计算量影响最大。

更大的步长（4 fs）可行，但需要用约束算法（LINCS/SHAKE）把所有含氢键长冻结——每步少算一些非键相互作用，代价是失去高频振动信息。团队在溶菌酶的折叠/去折叠模拟中比较过2 fs不冻结 vs 4 fs冻结氢：4 fs冻结氢的RMSD演化轨迹和2 fs的几乎重合，偏差<0.5 Å，但总模拟时间从14天缩短到7天。

制约步长上限的不是算法稳定性，而是氢原子的高频运动。如果把所有氢的质量重新分配到重原子上（HMR，Hydrogen Mass Repartitioning），振动频率降下来，步长可以推到4 fs甚至5 fs。不过HMR会略微改变动力学——对平衡态采样影响很小，但对非平衡动力学（如拉伸MD）要谨慎。

系综耦合：恒温恒压下怎么控

NPT系综是模拟生物大分子在水溶液中的标准方案。恒温控制常用Nosé-Hoover或Berendsen，Berendsen收敛快但不产生正确的正则系综分布，Nosé-Hoover是正确的但振荡收敛慢。团队的习惯做法是先跑500 ps的Berendsen做快速平衡，再切到Nosé-Hoover做生产运行——这样做兼顾了速度和正确性。

恒压控制是另一个容易出问题的地方。Parrinello-Rahman压浴加上各向同性耦合（isotropic），在膜蛋白模拟里可能会压偏——膜的法向（z）和膜平面（x-y）的压力分量本应独立控制。一个教训是膜蛋白OmpF的模拟：各向同性压浴下z方向的cell尺寸从117 Å跑偏到122 Å，换算成面积/脂质从62缩到55 Ų——这个偏差足以影响通道的开放状态。切到各向异性（semiisotropic，x-y耦合、z独立）后，z向cell尺寸稳定在实验值附近。

这套力场→时间步长→系综→恒温→恒压的耦合决策链，在全原子分子动力学中每一步都受前一步牵制。回到根本上，全原子分子动力学模拟原理的核心不是选对参数，是搞清楚每个参数决定什么——物理约束和数值约束之间的平衡，才是MD模拟真正的”原理”所在。体系特化参数的调试思路，站上有系列复盘可以查阅。

更多内容请访问 https://www.keyanxueshu.com/