第4章非键相互作用

4. 非键相互作用：范德华力、静电相互作用与截断处理

非键相互作用，说白了就是分子中那些没有直接成键的原子之间的「远程交流」。在分子动力学模拟里，这部分计算量最大，也最容易出问题。我个人习惯把非键作用分成两块来看：一个是范德华力，一个是静电相互作用。这两兄弟处理不好，你的模拟结果基本就废了。

范德华力，本质上是一种短程的、随距离变化很快的相互作用。它包含吸引和排斥两部分。吸引来自瞬时偶极的诱导效应，排斥则是因为电子云重叠时的泡利不相容原理。

描述范德华力最经典的模型就是Lennard-Jones势，简称LJ势。它的数学形式很简单：

V(r) = 4ε [ (σ/r)^12 - (σ/r)^6 ]

这里有几个关键参数：

你想想看，当r很大时，两项都趋近于零，没有相互作用。当r等于σ时，势能为零。当r小于σ时，r^(-12)项迅速占主导，产生强烈的排斥。当r在某个最佳距离时，势能达到最小值-ε。

核心要点：LJ势的r^(-12)项其实是为了计算方便而选用的，并非物理上的真实形式。真正的排斥项应该是指数形式的，但用12次方可以省去指数运算，计算效率更高。这是典型的「用精度换速度」的工程思维。

我在项目中遇到过一个问题：用OPLS力场模拟液态烷烃时，密度总是偏低。后来排查发现，是LJ参数的σ值偏大了0.02 nm。别小看这0.02，整个体系的密度直接掉了3%。所以参数选择一定要谨慎。

静电相互作用，说白了就是带电粒子之间的库仑力。它的形式比LJ势更简单：

V(r) = (1 / 4πε₀) * (qᵢ qⱼ) / r

其中qᵢ和qⱼ是原子所带的电荷，ε₀是真空介电常数。这里有个关键点：库仑势是长程相互作用，它随距离衰减的速度很慢（只有1/r）。

为什么会这样？因为LJ势的r^(-6)项衰减很快，到了1 nm以外基本可以忽略。但库仑势的1/r衰减，在2 nm、3 nm甚至更远的地方仍然有不可忽略的贡献。这就带来了一个很头疼的问题：截断半径怎么选？

注意：千万不要对静电相互作用直接使用简单的截断处理。我曾经见过有人把库仑势的截断半径设为1.2 nm，结果模拟出来的水溶液结构完全不对，介电常数都算错了。直接截断会在截断边界处产生人为的势能跳跃，导致能量不守恒。

截断半径，就是设定一个距离阈值，只计算这个距离以内的非键相互作用。对于LJ势，由于它衰减快，通常截断半径设在1.0-1.4 nm之间就足够了。但即使这样，截断后也需要做能量和力的校正。

对于LJ势的截断校正，常用的方法是：

对于静电相互作用，情况就复杂多了。常用的长程处理方法有：

我的建议：对于生物大分子体系，比如蛋白质、核酸，直接用PME方法。对于小分子液体或气体，反应场法可能就够用了。但不管用哪种方法，一定要做收敛性测试——改变截断半径和网格密度，看看能量和压力是否稳定。

下面这张图展示了非键相互作用的核心逻辑和参数选择流程：

说了这么多理论，最后给几条实战经验：

LJ参数组合规则要统一：不同力场用的组合规则不一样。CHARMM和AMBER用Lorentz-Berthelot规则（算术平均σ，几何平均ε），OPLS用几何平均。混用会导致参数不一致，模拟结果不可靠。
截断半径不要小于1.0 nm：我见过有人为了省计算量，把截断半径设到0.8 nm。结果模拟出来的液体结构完全不对，径向分布函数的第一峰位置都偏了。省下的计算时间全花在调试上了。
长程校正一定要做：对于LJ势，不做能量和力校正，模拟体系的压力会偏大。对于静电相互作用，不做长程处理，体系的结构和动力学性质都会失真。
验证三步走：先跑短时间的NVT模拟，检查能量是否稳定；再跑NPT模拟，检查密度和压力是否与实验值吻合；最后算径向分布函数，看局部结构是否正确。

避坑指南：我曾经在模拟离子液体时，直接用了水溶液的LJ参数。结果模拟出来的密度比实验值高了10%。后来才发现，离子液体的离子半径大，σ值需要重新拟合。所以，参数选择一定要基于目标体系的物理化学性质，不能盲目套用。

嗯，非键相互作用这部分内容就讲到这里。参数选择没有万能公式，关键是要理解每个参数的物理意义，然后根据你的体系特点做针对性调整。多跑几个测试，多对比实验数据，慢慢就能找到感觉。