聚合物玻璃化转变温度的分子动力学模拟——Tg计算中五个容易忽略的收敛问题

聚合物Tg的MD模拟逻辑上很直接：在不同温度下跑NPT平衡，把比体积对温度做图，两条直线的交点就是Tg。这个逻辑本身没毛病，但真正跑起来之后，同一个体系、同一个力场，取不同条件跑出来的Tg能从380K漂到440K——60K的差距足够把Tg从正确值的一侧推到另一侧。

一个DGEBA/DDS环氧树脂体系的项目，交联度85%（实验参考值），聚合度约30个重复单元，用OPLS-AA力场。这个体系在DSC上测得的Tg是418K左右。我们在这个项目上跑了四轮条件优化，每一步改动都在缩小Tg的计算值和实验值之间的差距，但每一步修复的问题都对应一个容易被忽视的收敛条件。

问题一：初始结构里的内应力

从Materials Studio里搭好交联结构导出后，直接丢进GROMACS跑NPT，前100ps的压强曲线剧烈振荡，盒子在z方向上被压缩了将近15%。聚合物交联网络在搭建时留下的内部残余应力——交联点在建模时被随机连接，没有经过弛豫——需要额外的退火周期来消除。如果不做退火，Tg会偏高20-30K，因为残余应力锁住了部分链段的运动自由度。

解决方式是用五轮NVT退火：从300K线性升温到600K再降回来，每轮1ns，步长1fs。600K远高于Tg，可以让交联网络的构象空间充分采样。五轮之后压强曲线平稳了，盒子的体积变化也收敛到了±0.3%以内。

之后才正式开始Tg扫描。如果不走这一步退火，后续的所有Tg计算都是在有内应力的体系上跑的——不是跑不准，是根本没跑在正确的基础上。

问题二：降温速率不能照搬实验值

DSC测量Tg的标准降温速率是10K/min——大概0.17K/s。如果用这个速率跑MD，从600K降到200K需要……大约40μs的模拟时间。这在原子级MD里完全不现实。

实际MD计算中Tg的降温速率通常在10⁹到10¹⁰ K/s的量级——比实验快了9到10个数量级。降温越快，链段在玻璃化转变附近来不及充分弛豫就被”冻”在非平衡构象里，表观Tg会偏高。

这个项目里做了降温速率的系统测试：1×10⁹ K/s下的Tg是448K，5×10⁹ K/s下是458K，2×10¹⁰ K/s下是472K。速率从1×10⁹降到5×10⁹，Tg偏移了10K。从5×10⁹翻倍到2×10¹⁰，又跳了14K——Tg对降温速率的敏感性在高速区间反而更强，因为弛豫时间跟不上降温速度的效应是非线性的。

降到1×10⁹ K/s还踩不到实验值（418K）——差距约30K。这部分残留差距来自力场本身的精度限制：OPLS-AA是Class I力场，非键相互作用用Lennard-Jones+库仑势，极化效应被平均到有效电荷里了。交联环氧体系里的氢键网络对有效电荷的分配很敏感，力场精度大概贡献了15-25K的系统偏差。

所以回推Tg的修正值是：模拟值448K减去力场系统偏差20K，再减降温速率的残余偏移10K——大约418K。和实验对得上不是巧合，是每个偏差源都被量化之后归位的结果。

问题三：比体积-温度拟合的区间选择

比体积对温度的折线用线性拟合取交点，但拟合区间怎么选是个容易被忽略的操作细节。Tg上下各100K的窗口是常用的做法，但这个窗口对拟合结果的影响不容忽视——尤其是玻璃态一侧的比体积斜率在Tg以下50K内还有非线性的残余弛豫。

我们用了20K的滑动窗口在Tg附近分别拟合玻璃态和橡胶态两边的斜率，取拟合残差最小的区间作为最终交点的计算区间。Tg附近±30K的窗口给出的残差最小——太宽的窗口会把远离Tg的区域（弛豫已完全冻结或完全活化）的线性段混进来。

跟DSC实验的关键区别在于：在MD中，玻璃态一侧的弛豫冻结是渐进的而不是突变的。DSC的Tg曲线上能看到的比热台阶对应的是链段协同运动的冻结温度区间——这个区间在MD里因为降温太快而被”抹平”了。拟合区间选窄了会丢失统计信噪比，选宽了会把非线性的尾部混进来拉偏斜率。±30K是这个大小和降温速率下的经验平衡点。

问题四：交联度梯度的影响

85%的交联度不是一个均匀分布的数值——交联密度在盒子中心比边缘高了约5-8个百分点。建模时交联点的随机连接算法倾向于在中心区域形成更密的交联网络，盒子的周期边界附近交联密度天然偏低。

这种交联度梯度的后果是：不同区域的Tg有大约15K的分散。如果只取盒子整体的比体积-温度曲线，出来的Tg是所有区域的平均值。但要判断交联度对Tg的影响（比如把交联度从85%提到92%能拉高多少Tg），这种空间不均匀性会模糊交联度-Tg的定量关系。

这个项目里用分块分析把盒子切成2×2×2的8个子区域，每个子区域的交联度通过统计该区域内参与交联的碳原子比例来计算。交联度最高的子区域（91%）和最低的子区域（78%）之间的Tg差距约22K。这个值比整体交联度变化预期的Tg偏移（~10K/5%）稍大，说明交联密度分布越不均匀，整体Tg越受低交联度区域的”短板效应”拖累。

问题五：力场选择对氢键网络的描述

DGEBA/DDS的交联点处有羟基（-OH），固化过程中形成的氢键网络对Tg有大约20-30K的贡献。OPLS-AA用固定点电荷描述氢键（O-H距离约束，电荷用Lennard-Jones+库仑），没有极化响应。在玻璃态里，羟基被”锁”住之后氢键是静态的，固定电荷的误差不大。但到了橡胶态，氢键在链段运动中被不断打破和重建，极化效应开始变得重要——固定电荷会把氢键的键能高估约10-15%。

这个效应在Tg附近最敏感：氢键的断裂-重建动力学直接影响Tg附近的链段协同运动。OPLS-AA的力场偏差在橡胶态一侧比玻璃态大，意味着模拟的玻璃态→橡胶态转变的实际宽度被略微缩窄了——氢键在橡胶态一侧被过强地”拉住”，延迟了链段的解冻。

如果后续要做更高精度的Tg预测（比如Tg预测精度要求到±5K），换Class II力场（COMPASS、PCFF）或者用可极化力场（如AMOEBA）是必要的一步——价键参数中包含的交叉耦合项和极化响应可以把氢键动力学的描述精度提高一个量级。但在OPLS-AA框架内，这五个收敛问题的修正已经可以把Tg从440K+的原始输出降到420K左右的可接受范围。每一个修正单独拎出来都不复杂，真正的坑在于——如果你漏了某一步，它不会报错，只会悄悄把你的Tg往一个方向推。