第三章 高分子力场基础

做高分子模拟这些年，我最大的体会就是：力场选对了，模拟就成功了一半。力场这东西，说白了就是描述原子间相互作用的数学函数。你给它一套坐标，它告诉你原子间有多大的力、多少能量。嗯，就这么简单。

但真正做起来，你会发现力场里的门道特别多。今天我就把高分子力场里最核心的几个部分拆开来讲：键伸缩、键角弯曲、二面角扭转，还有非键相互作用。这些都是你搭建高分子模型时必须面对的基础问题。

3.1 键伸缩能

先说说最简单的——键伸缩。两个原子通过化学键连接，就像用弹簧连着一样。你拉长它，它想缩回去；你压短它，它想弹回来。这个能量变化，我们用谐振子模型来描述：

E_bond = ½ × k_b × (r - r₀)²

其中 k_b 是力常数，r₀ 是平衡键长。我在做聚乙二醇模拟时遇到过一个问题：默认的平衡键长参数跟实际结构对不上，结果链构象完全跑偏了。后来我学乖了，每次建模前都会先查一下目标分子的键长数据。

这里有个坑要提醒你：

3.2 键角弯曲能

三个原子形成一个键角，比如C-C-C。这个角度的弯曲也需要能量。公式跟键伸缩类似：

E_angle = ½ × k_a × (θ - θ₀)²

θ₀是平衡键角，k_a是角度力常数。你想想看，如果键角参数给错了，整个链的局部构象都会变形。我做聚苯乙烯模拟时，苯环上的C-C-C键角参数稍微偏了一点，结果侧基的取向全乱了。

我个人习惯的做法是：先用小分子模型做气相优化，把键角分布跟实验或高精度量化计算结果对比一下。没问题了再上高分子体系。

3.3 二面角扭转能

这个就有点意思了。四个连续原子形成的二面角，决定了高分子链的局部构象。比如聚乙烯的trans和gauche构象，就是由二面角决定的。

二面角能量用傅里叶级数表示：

E_dihedral = Σ [V_n × (1 + cos(nφ - γ))]

V_n是势垒高度，n是周期性，γ是相位角。为什么用傅里叶级数？因为二面角势能面通常有多个极小值点，一个余弦项搞不定。

我记得有一次做聚丙烯模拟，默认的二面角参数导致gauche构象占比过高。查了半天才发现，是参数里少了一个n=3的项。加上之后，构象分布就跟实验对上了。

3.4 非键相互作用

非键相互作用包括两部分：范德华力和静电作用。它们描述的是没有直接化学键连接的原子之间的作用。

3.4.1 范德华力（LJ势）

最常用的是Lennard-Jones 12-6势：

E_LJ = 4ε × [(σ/r)¹² - (σ/r)⁶]

ε是势阱深度，σ是碰撞直径。r⁻¹²项描述短程排斥，r⁻⁶项描述长程吸引。说白了，两个原子靠太近会排斥，离太远又没作用，中间有个最佳距离。

做高分子熔体模拟时，LJ参数的截断半径怎么选很关键。我一般用1.2 nm，但有些体系需要更大。你想想看，如果截断半径太小，长程的色散作用被忽略了，密度和压力都会算不准。

3.4.2 静电作用（库仑势）

E_coul = (q_i × q_j) / (4πε₀ × r)

q是原子电荷，ε₀是真空介电常数。静电作用在极性高分子中特别重要，比如聚酰胺、聚酯这些。

电荷怎么分配？这是个老大难问题。我试过Mulliken、NPA、RESP好几种方法，最后发现RESP电荷在凝聚相模拟中表现最好。但要注意，RESP电荷计算成本高，大分子体系得用分段拟合。

3.5 力场参数化方法

力场参数从哪来？主要有三条路：

1. 从已有力场移植：比如OPLS、CHARMM、AMBER这些通用力场，参数都是现成的。适合常见高分子。
2. 量子化学计算拟合：用小分子模型做DFT计算，然后拟合力场参数。精度高，但费时间。
3. 实验数据反推：用晶体结构、振动光谱、热力学数据来优化参数。适合特定体系。

我个人最常用的是第二种。先做DFT优化，得到平衡几何和势能面，然后用遗传算法或模拟退火来拟合参数。这里有个经验：拟合时不要只看能量，还要看力和Hessian矩阵，这样参数更鲁棒。

下面这张图总结了高分子力场的核心组成和参数化流程：

最后说一句，力场参数化是个细致活。我见过不少新手，随便找个力场就开始跑模拟，结果算出来的东西跟实验对不上，还以为是自己的方法有问题。其实很多时候，问题就出在力场参数上。

做模拟之前，花点时间验证力场参数，这笔时间绝对值得花。

---