目录

使用 APPLE&P 计算离子液体的粘度

Apple&P是一个极化力场,可用于带电体系,如电解质(例如电池中的电荷迁移率)、离子液体类体系的分子动力学模拟。

在本教程中,我们将以 N-甲基-N-丁基吡咯烷鎓 (pyr14) 和双(三氟甲磺酰亚胺) (TFSI)(1:1 比例)的混合物为例,描述如何使用 APPLE&P 极化力场计算电解质混合物的粘度。

版本需求: AMS2023 或更高版本。

注意:

建模

下载并解压,得到pyr14、TFSI分子的xyz文件。

创建混合物模型:

注意:此时的密度很低,与实际密度不一致。这是很常用的一种方式,因为体系要求分子保持一定间距,而这种建模方法,分子是刚性堆叠的,因此不一定能堆叠出高密度的体系。因此一般先创建低密度体系,然后逐渐压缩到特定的密度。

如果使用 Python 脚本堆叠建模代码,参考:Python 堆叠建模

第一步:预优化、压缩密度

切换到 ForceField 模块,并设置参数如下: 点击Type: Apple & P后面的 > 按钮,为两种离子设定电荷,并点击 Generate 按钮: 此时,APPLE&P 工具将进行连通性分析并识别系统中每个元素的原子类型,然后在您的工作目录中写入一个扩展名为.ff的力场文件。如果输入电荷不正确,AMSinput 窗口将报错,成功后,窗口左下角出现:Finished generated APPLE&Patom types 的提示内容。

这部分工作,对应Python脚本中的(如果不使用Python则可跳过不看):

s = Settings()
s.input.ForceField.Type = 'APPLE&P'
s.input.ForceField.ForceFieldFile = os.path.abspath(forcefield_file)

预优化

因为密度尚未正确,因此这里选择优化Lattice,Convergence设置为Basic大致优化一下就可以了。

保存作业并运行。运行完毕,图形窗口提醒是否更新为最新优化后的结构,点Yes即可。

压缩

将Main窗口的Task 改为Molecular Dynamics:

设置NVT系综: 这里Checkpoint频率设为10000,是避免生成大量中间文件。一般来说都可以把这个值设置的很大。起始温度333K,如果不设置将沿用恒温器设置里面出现的第一个温度,因此不设置也可以(333与300K的差异并不影响最终结果)。

然后设置压缩:将Target length设置为 25Å、25Å、25Å,即如此最终密度将接近实验密度(用户在Builder中建模时,分子填充完毕后,可以修改Lattice vector看看晶格矢量为多大的时候,密度是合适的,记住那个值,然后改回原值35×35×35)

再次保存覆盖之前作业,并运行之。在SCM → Movie 中可以看到体积的变化,Movie → MD Properties → Density可以看到密度的变化。运行完毕后,体系结构更新,待下一步计算。

第二步:NVT系综热化平衡

首先将体系在NVT系综下进行热化达到平衡。分子动力学参数设置可以沿用之前的设置,但需要点击Model → MD Deformations中Deformation:#1参数块前面的 - 以去除形变参数。这里仍然运行了10000步,但您的体系可能需要更长时间,能量上逐渐达到平衡为最佳。

保存作业并运行。运行完毕后,图形窗口仍然弹出窗口寻味是否更新结构,建议同时勾选 Use MD velocities from result 后点击Yes。然后进行下一步。

第三步:NPT系综下平衡模拟

在原先设置的基础上,增加恒压器设置(总步数可以设置的较大,以便体系达到平衡):

另存作业或覆盖前面作业,并运行。

单击SCM → Movie来分析此次 NPT 作业。Movie → Graph → Add Graph,然后MD Properties → Density观察密度的变化。理想情况下,NPT 收敛后的密度与实验较为接近,这里我们发现 NPT 模拟得到的密度在 1336 至 1346 kg/m$^3$之间,比实验密度 1367 kg/m$^3$小约 2%。

运行完毕后,图形窗口仍然弹出窗口寻味是否更新结构,建议同时勾选 Use MD velocities from result 后点击Yes,然后进行下一步。这里,如果体系足够大,能观察到压强的震荡,对应的也有密度的震荡,在已经达到平衡的区域,也可以找到压强最准确的一帧,更新结构进行后续的计算。

第四步:正式计算粘度的分子动力学模拟

去掉前一步的压强设置(Barastat改为None),步数改到非常大的值,例如400000,Sample frequency可以改到较大例如1000,Checkpoint frequency改为400000。一般建议模拟时长到几ns,因此远大于此处演示所采用的40万步。

点击Trajectory Options,勾选Calculate pressure:

Details → Expert AMS 下拉到 MolecularDynamics Binlog 勾选 Time 与 Pressure tensor(将在模拟的每一步写入压力张量到输出文件)。

保存或另存作业并运行之。

粘度的结果分析

电解质的粘度可以使用 Green-Kubo 关系式计算:

\[\eta = \frac{V}{k_{\rm{B}}T} \int_0^\infty \langle P_{\alpha\beta}(0)P_{\alpha\beta}(t) \rangle \rm{d} t\]

其中 η 是粘度,V 是细胞体积, k$_B$ 是玻尔兹曼常数,T 是温度, P$_{αβ}$ 是压力张量的非对角分量,t 是时间,\(\langle P_{\alpha\beta}(0)P_{\alpha\beta}(t) \rangle\) 是一个自相关函数,其平均值取自所有时间原点和非对角线分量(αβ为yz、xz或xy)。

用户可以使用以下 Python 脚本(下载)后解压,从ams.rkf文件中的 binlog 段读取压力张量,并根据非对角压力张量自相关函数的积分计算粘度。所得粘度值为双指数拟合的极限值A ,定义为:

\[f(x) = A\left( \lambda (1-\exp(-x/\tau_1))) + (1-\lambda)(1-\exp(-x/\tau_2))\right)\]

脚本运行方式:下载后解压到作业所在目录,与*.results文件夹同级。AMSJobs也进入该目录,Help → Commandline输入sh回车,然后输入命令:

$AMSBIN/amspython viscosity_extractor.py *.results/ams.rkf回车

例如:

注意

如果最后一步的分子动力学模拟长度为几纳秒,则 Python 脚本的执行可能需要长达几个小时。事实上,对数百万帧的自相关函数进行评估非常耗时。

得到:

基于此处 400 ps 动力学过程的模拟,我们得到的极限值为 η= 23.92 mPa·s。请注意,尽管此模拟很短,但我们与混合物的实验粘度值 21.0 mPa·s(在 333 K 下实验值,参见O. Borodin, Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids, J. Phys. Chem. B, 2009, 113, 11463–11478)非常吻合。不过,粘度的波动表明模拟没有很好地收敛。

如果我们继续之前的模拟直至 5 ns,我们会得到一个极限值 η= 22.55 mPa·s。在这里我们可以很容易地看到粘度收敛,趋向稳定值:

总结与经验