目录

使用 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。在这里我们可以很容易地看到粘度收敛，趋向稳定值：

总结与经验