在本教程中,我们将利用 ForceField 引擎进行分子动力学模拟,计算水溶液中氯化钠离子的离子电导率。
这里介绍的体系和 workflow 与文献中描述的非常相似:
本教程包括以下步骤:
AMSinput → ForceField模块 → Type: → Amber95。创建一个Na原子,点击Atom Type右侧的>按钮,对Na原子进行Amber原子类型的设置,Na类型为IP,电荷为1:
File → Export coordinates → .in,保存为Na.in。
类似地创建Cl离子,FF Type为IM,电荷为-1,导出为Cl.in;创建水分子:FF Types O: OW, H: HW, FF Charges O: -0.8340, H: 0.4170,并导出结构为h2o.in。
在菜单栏Edit → Builder中创建这三个组分的均匀混合物:
注意调整晶格尺寸,让New density will be around显示密度为1.029g/ml左右,如此会产生密度为 1.029 g/mL 的 NaCl 水溶液体系类似上图所示。
注意创建体系时,请先确保力场选择的是Amber95。
首先粗略优化体系:
Main → Task → Geometry Optimization,确保Periodicity为Bulk,保存作业并运行。运行完毕后,询问if the structure in AMSinput should be updated时,点yes,将优化的结构更新到AMSinput窗口。
然后进行弛豫:
运行结束,SCM → Movie → Properties → Potential Energy,再Properties → Temperature,
可以看到,体系的势能在模拟过程中,逐渐收敛到一个稳定值,温度在 300K 上下振荡。同样地将最后的结构更新到AMSinput中。
Task → Molecular Dynamics,点击Task → Molecular Dynamics右侧的>按钮,进入MD详细参数设置。
点击Thermostat右侧>按钮设置恒温器:
File → Save as另存一个作业,例如NaCl_production,并运行作业。
离子电导率与原子电荷加权均方位移的斜率成正比:
\(\sigma = \lim_{t\to\infty} \frac{1}{2tdVk_BT} \sum_i^N q_i^2 \langle ({\bf{r_i}}(0) - {\bf{r_i}}(t))^2 \rangle\)
这里d是系统的维数(通常d=3), N是离子的数量,Movie中可以计算σ值及其斜率(电导率)。
AMSMovie → MD Properties → MSD:
点击Generate MSD按钮,Movie窗口中出现两条新曲线。底部曲线描绘了离子电导率随时间的收敛情况,离子电导率的最终值为7.34 S/m(或 C²/Jms)。文献指出,该值应介于 1.2 和 13.5 S/m 之间,多次模拟的平均值是 5.9 S/m。
顶部曲线描绘了电荷加权均方位移曲线,离子电导率是其斜率。该图显示,在较短的相关时间下,曲线的斜率很大,然后逐渐变为恒定。曲线变为线性的点大致对应于离子运动“去相关”所需的时间,Movie会根据均方位移曲线的形状来估计这一点。实际上,底部曲线是从大约 7ps 开始的,这是基于 Movies 的这种估计。
不过,肉眼简单看一下顶部曲线,可以看到,“去相关”时间的更好猜测是在 10 ps 左右。我们可以根据这个新的“去相关”时间重新计算离子电导率。
MD Properties → MSD 保持之前的设置,然后增加设置Start time slope为10000.0(注意单位为fs):
点击Replace MSD,离子电导率的新值为7.29 S/m。