LAMMPS计算结合能：聚合物-纳米填料界面的结合能，从拔出模拟到PMF，力场精度决定你拉出来的是多少

LAMMPS计算结合能——尤其是聚合物-纳米填料界面的结合能——定义上很清楚：把一个纳米填料从聚合物基体中拉出来所需的功除以界面面积。但怎么在线性MD模拟中”拉出来”，有三种操作方法，算出来的数字可以差一倍。

一个做CNT增强环氧树脂的项目，碳纳米管直径约1.2nm（(10,10)扶手椅型），埋在交联的DGEBA/DETA环氧基体中。实验上用AFM的pull-out测量，单根CNT的界面剪切强度约40-60MPa。项目要做的就是用MD复现这个数量级，然后筛选不同的CNT表面官能团化对结合能的增强效果。

方法一：静态减法

最直接的方法——算一个总能量E_total（聚合物+CNT体系）、一个聚合物的能量E_polymer（去掉CNT但保留空洞）、一个CNT的能量E_CNT，结合能 E_binding = (E_polymer + E_CNT − E_total)。

但聚合物在CNT拔除后的空洞不会保持原样——在常温MD里空洞周围的聚合物分子会立刻重排、坍塌甚至部分填充空洞。E_polymer取的是带空洞的构型的能量，而这个构型在有限温度下不稳定。静态减法的E_binding在0K下有意义，在300K下没有热力学定义——结合能不是焓变，是自由能变。

这个项目里静态减法算出来的CNT/环氧结合能是2.3 J/m²，对应界面强度约160MPa——是实验的3-4倍。不是CNT和环氧真的结合那么强，是聚合物空洞的应变能没有被从总能里正确排除。

方法二：受控拔出模拟

升级方案是直接模拟CNT的拔出过程：固定CNT的一端，以恒定速度（0.1 Å/ps）沿轴向把CNT从环氧基体中拉出来。记录拔出过程中的轴向力-位移曲线，积分曲线下的面积就是拔出功。

力-位移曲线上有一段典型的尖峰（初始脱粘）——最大拔出力的存在来自CNT和聚合物之间的范德华和静电力累积。对于无官能团化的CNT，最大拔出力在LAMMPS里是约0.92 nN。对CNT-COOH（羧基官能化，覆盖率约5%），最大拔出力升到1.41 nN。

但拔出速度是关键变量。0.1 Å/ps这个速度在物理上快得离谱——相当于CNT以10m/s的速度被拉出来。在这种速度下，聚合物链没有时间松弛——拔出力被”粘性阻力”的成分污染了。将速度降低到0.01 Å/ps和0.001 Å/ps，最大拔出力分别降到了0.78 nN和0.65 nN——接近速度无关的平台（准静态极限）。

0.1 Å/ps的”快速拔出”高估了结合能约40%。

方法三：伞形采样构建PMF

从准静态拔出到热力学上最严格的方法，是用伞形采样（Umbrella Sampling）沿拔出反应坐标ξ（CNT质心相对于聚合物基体界面的位移）构建平均力势（PMF）。

在ξ从0到20Å的范围内设置40个窗口，每个窗口用谐波约束（力常数3 kcal/mol·Å²）限制CNT在±0.25Å范围内波动。每个窗口采样5ns，40个窗口×5ns，共200ns。用WHAM算法重构PMF剖面。

PMF的拔出力势垒高度是1.85 eV/nm（对应约0.30 nN拔出力）——比准静态拔出模拟的0.65 nN还低了一半。为什么？

因为伞形采样考虑了聚合物链在CNT拔出过程中的构型重排——聚合物链在ξ的每个位置上都是在局部热力学平衡的。准静态拔出虽然速度低了，但每个位移位置没有给足够的构型采样时间——聚合物链的部分构型被”冻结”在了拔出路径上。

PMF给出的界面剥离能是0.72 J/m²，对应界面剪切强度约48MPa——落在AFM pull-out实验40-60MPa的范围内。

精度和代价

从静态减法→拔出模拟→PMF，结合能从2.3→1.2→0.72 J/m²，每次”升级”都降了一大截。这不是方法在变，是对”结合能”这个物理量的定义在越来越接近真实——从焓变到自由能变。

代价也是递增的：静态减法跑3个单点能（几分钟）；准静态拔出跑约20ns（几小时）；PMF完整剖面跑200ns（两三天）。选择哪种方法不取决于”哪个更准”，取决于你的科学问题需要自由能精度到什么程度。如果只是用结合能做定性排序（官能团A比B强多少），拔出模拟的速度依赖趋势已经够用了。如果要和实验的界面剪切强度做对标，PMF是唯一正确的答案。