使用 APPLE&P 计算离子液体的粘度
Apple&P是一个极化力场，可用于带电体系，如电解质（例如电池中的电荷迁移率）、离子液体类体系的分子动力学模拟。
在本教程中，我们将以 N-甲基-N-丁基吡咯烷鎓 (pyr14) 和双(三氟甲磺酰亚胺) (TFSI)（1:1 比例）的混合物为例，描述如何使用 APPLE&P 极化力场计算电解质混合物的粘度。
版本需求： AMS2023 或更高版本。
注意：
- 为了提取可靠的粘度值，体系必须达到很好的平衡状态。本教程仅仅为了演示功能，因此平衡时间非常短，实际计算建议在集群进行，令体系达到充分平衡
- 本教程除了通过下面的图形界面操作，也可以通过Python脚本一次性完成，脚本下载（点击）后解压即可直接通过命令行运行：amspython applenp_viscosity_md.py回车。参考：命令行运行 Python
建模
下载并解压，得到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。在这里我们可以很容易地看到粘度收敛，趋向稳定值：
总结与经验
- 模拟需要很好地平衡：本例中，我们从低密度开始，对模拟盒进行变形以接近实验密度，用户还可以考虑对系统进行温和的退火处理（关于退火的概念与操作，参考Li-S中Li离子的扩散系数：MSD法、VACF法，以及结合尔伦尼乌斯公式将扩散系数外推到低温中《退火》小节，该教程中只做了一次升温降温，实际上一般的退火会升温降温10做来个循环）达到更好的平衡。
- 一般来说，使用 NPT 模拟平衡液体的时候，要注意验证密度、压力是否接近实验值。
- 粘度通常对体系的规模大小很敏感，应使用更大的体系以获得更稳定的结果。
- 在整个时间范围内，粘度曲线应平滑上升，平滑特征代表收敛、数据的稳定性。
- 如果粘度曲线有噪声，请减少最大相关时间（Python脚本中的max_dt_fs，单位为fs）或运行更长时间的 MD 模拟。拟合的双指数曲线在任何点都不应大幅偏离粘度曲线。
- 确保模拟时间至少比最大相关时间长 20 倍（但最好接近至少 50 倍）。例如：如果 MD 模拟长度为 2 ns，则将后分析脚本中的 max_dt_fs 设置为 (2 ns / 20) = 100 ps 或更短。在上图中，最大相关时间太长（400 ps / 100 ps = 4 ），大量噪声也从侧面印证了这一点。
- 对积分的最大贡献来自较短的相关时间。因此，有必要使用 BinLog 功能，该功能在每个时间步写入压力张量。