热固性聚合物的玻璃化转变温度

在本ReaxFF教程中，解释了从聚合物密度的温度依赖性计算玻璃化转变温度T_g。 详细系统性研究，参考文献：

• Effect of chemical structure on thermo-mechanical properties of epoxy polymers: Comparison of accelerated ReaxFF simulations and experiments, Polymer 159, 354-368 (2018).

本教程包含如下几个步骤：

1. 导入聚合物结构
2. 模拟退火
3. 提取密度-温度性质
4. 计算玻璃化转变温度

本计算对硬件要求较高，建议在集群或工作站上运行，软件版本要求AMS2021以上。本文以ReaxFF作为范例，其他模块例如机器学习势、DFTB、基于BAND模块则属于AIMD，除了Main的参数设置不同，其他设置完全一致。

这里我们使用bond boost法生成的交联环氧聚合物，当然如果有其他结构的*.xyz文件，也可以直接通过AMSinput → File → Import Coordinates导入。

为了精确模拟密度的温度依赖性，我们必须确保聚合物结构尽可能均匀。聚合物结构内部的大孔是不均匀结构的指标，这将导致非常嘈杂和难以分析的结果。模拟退火是一种简单的优化策略，旨在在复杂的高维势能面上产生全局最小能量结构。在退火过程中，MD模拟过程中的温度逐渐升高，希望系统能够克服越来越多的能垒，有效探索能量势能面全局。加热阶段之后是冷却阶段，在此期间，系统被冷却回初始温度。通常，结果结构的能量明显低于起始结构，尽管不能保证它是最佳结构。

参数设置如下：

退火过程中，体系的体积有可能扩大或者缩小，因此使用恒压器（此处使用恒压器还有一个原因是原子的个数足够多，如果体系只有几百元子、一两千原子，恒压器效果就不那么好，原子个数越多，恒压器控压的效果越好）

接下来，我们定义所需的温度程序。为了从计算中提取密度，我们将把温度从298.15提高到598.15，这对于模拟退火来说仍然很低，但对于我们的目的来说已经足够了。我们要实施的温度程序很简单：花费30000步升温25K，随后10000步恒温对密度进行取样，循环此步骤直到达到598.15 K，然后反向冷却至298.15 K。如下图所示意（下图仅仅示意了两步升温、采样，实际需要非常多步升温、采样，详见下文中温度的设置即可知）：

在“Thermostat”中设置温度：

298.15 298.15 323.15 323.15 348.15 348.15 373.15 373.15 398.15 398.15 423.15 423.15 448.15 448.15 473.15 473.15 498.15 498.15 523.15 523.15 548.15 548.15 573.15 573.15 598.15 598.15 573.15 573.15 548.15 548.15 523.15 523.15 498.15 498.15 473.15 473.15 448.15 448.15 423.15 423.15 398.15 398.15 373.15 373.15 348.15 348.15 323.15 323.15 298.15 298.15

“Duration”中设置每个阶段的步数：

10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000 30000 10000

如下图所示：

保存作业并运行，本作业约8000原子，建议使用8~16核计算。

为了从轨迹中提取密度和温度进行后处理，我们使用Python脚本densities.py。将其与*.ams输入文件放在同一文件夹中。

作业运行完毕后，运行该Python脚本提取密度-温度数据：

如果使用Windows系统：

1. 双击AMS202*.*/ams_command_line.bat，弹出命令行窗口
2. 输入bash回车
3. 进入作业路径，例如cd D:/ADF_DATA/AMS_Task/SimulatedAnnealing/回车，注意斜杠方向和Windows显示文件夹路径的斜杠方向相反
4. 输入plams densities.py -v resultsdir=*.results回车，*号替换为作业名字

如果使用Linux系统：

1. 进入作业路径，例如cd /home/user01/SimulatedAnnealing/回车
2. 直接输入plams densities.py -v resultsdir=*.results，*号替换为作业名字

脚本运行后，会在该窗口生成五列数据，分别是：温度、密度、晶格常数a、b、c。

将该温度、密度数据使用作图软件，例如Excel、Origin做出（如下左图所示）：

包含两条曲线，分别对应升温和降温过程。

如果初始结构和平衡结构的最终密度显著不同，且实验密度（如果可用）相差很远，则密度不会收敛，应进行另一次模拟退火。如果最终密度和初始密度接近，则可以认为模拟是收敛的(如上右图所示)，可以继续计算T_g。

为了计算Tg，通过下图所示的数据点进行两次线性拟合。这两条线之间的交点标记玻璃化转变温度。

应选择适合线性拟合的子集点，以便获得最大值（r^2）。如果r^2值太低(即数据太嘈杂)，或曲线不相交，我们建议尝试以下故障排除方法之一：

• 运行另一个模拟退火模拟
• 扩大温度范围
• 减少或增加densitions.py脚本中的平均参数

本文参考英文教程：https://www.scm.com/doc/Tutorials/MolecularDynamicsAndMonteCarlo/PolymersGlassTransitionTemp.html