因为粘度的计算不涉及反应,因此除了使用ReaxFF外,也可以使用Force Field模块中其他力场(例如GAFF力场)或机器学习势,甚至也可以使用第一性原理进行分子动力学模拟,从而类似地生成粘度信息。分子动力学的参数设置大同小异,只是Main面板的参数设置有所差别。
这里以苯为例,室温下苯的粘度实验值为0.5 mPa s。本文计算粘度是通过分子动力学,因此单组份、多组分混合物都是适用的,文末介绍了使用QSPR方法计算粘度,只适用于单组份液体。
AMSinput - Edit - Builder,设置晶格常数为26.666*26.666*26.666 Å3,并加入128个Benzene分子(创建分子混合物参考教程:【入门基础教程】创建均匀混合物),从而密度与实验密度尽量接近:
注意体系的大小,会严重影响粘度的计算数值,因此体系规模应该大一些,这样粘度数值可靠性强一些(统涨洛误差减弱)。
基本的设置(周期性一定要确保为Bulk): 由于不关心化学反应,因此步长不需要太小,1fs就可以了,但是运动轨迹对粘度非常重要,因此每隔5fs就需要保存一次轨迹(甚至可以更小),关于动量、角动量守恒的限制也可以去掉(起始温度也可以不设置,因为后面Thermostat会设置温度): 温度的设置: 这里不建议使用NPT系综(因此不需要设置Barostat),因为Green-Kubo公式中有体积,如果使用NPT系综,则体积会发生变化,从而导致求解非常麻烦,另外液体在特定温度下,密度是确定的,所以使用NVT系综即可。
保存作业并运行。在logfile里面,确保Pressure一列有具体数值。如果使用ForceField模块,则未必有这个数值,有可能是“n/a”,那就需要增加设置:Model → MD → Trajectory Options后面的详细设置按钮 → 勾选Calculate pressure(为了减小计算量,可以去掉Write velocities、Write charges、Write bonds、Write molecules选项)。
压力张量自相关函数的计算范围,必须是动力学模拟达到平衡之后的时间。所以如果体系事先经过弛豫,则除去前面少量帧数,后面大量帧数都可以使用,如果没有弛豫,在势能达到稳定前的帧数,是不能用的。
从Movie中显示的势能的结果来看,300帧的时候,就接近平衡状态了,因此我们可以从300帧开始生成自相关函数(总共12001帧,模拟时间太短,帧数太少,会导致粘度数据不可靠):
MD Properties - Autocorrelation Function…
点击Generate ACF按钮,即可生成相关的几个函数,这里我们需要关心的是其中的积分函数,这个函数就是Green-Kubo公式中的积分。如下图绿色框所示,积分值为1.25283 × 10-10 Hartree2fs/Borh6:
带入Green-Kubo公式,并统一到国际标准单位:
因此粘度η = [18872 × 10-30 / (1.38 × 10−23 × 298) ] × 1.25283 × 10-10 × (2.942 × 1013)2 × 10-15 = 0.0005 Pa s= 0.5 mPa s
与实验符合的极好。
AMS2024版以后,该输出单位改为了Js/m3。
AMS的COSMO-RS模块,也能预测液体粘度:使用QSPR预测熔点、液体粘度、蒸汽压,以及三相点温度、合成难易程度、溶度参数、等张比容、介电常数等