自扩散系数计算

本文主要简介了通过GROMACS软件计算SPC/E水模型自扩散系数的步骤。

体系构建和top生成

通过GROMACS软件内置的gmx solvate可以很容易的构建水盒子,然后通过gmx pdb2gmx去产生SPC/E水模型的top文件。在这里,我采用的是OPLS-AA力场,构建了一个(4 x 4 x 4) $nm^3$的水盒子。完整命令依次如下:

1
2
gmx solvate -cs spc216.gro -box 4 4 4 -o water.gro
gmx pdb2gmx -f water.gro -o spce.gro -ff oplsaa -water spce

初始体系如下图:

模拟过程

能量极小化

首先通过能量极小化步骤,调整体系不合理的结构布局。

1
2
gmx grompp -f em.mdp -c spce.gro -o em.tpr
gmx mdrun -deffnm em -v

NVT模拟

经过能量极小化以后的体系就可以在NVT系综下进行模拟。因为我是采用的gmx solvate命令构建的水盒子,理论上此时的水密度已经是最合适的了,所以不需要再NPT系综下进行控压和体积的调整。如果你是采用的Packmol等软件通过指定水分子数目来构建的体系,那么此时的水密度显然不是最好的,因此建议先在NPT下进行体积调整,然后NVT模拟。
至于为什么要采用NVT进行成品模拟,主要是因为NPT系综压力对系统动力学的干扰。具体推荐你详读该文章Best Practices for Computing Transport Properties 1. Self-Diffusivity and Viscosity from Equilibrium Molecular Dynamics。了解了这些,你就可以接着做NVT模拟了,我这里模拟的温度为298K。完整命令如下:

1
2
gmx grompp -f nvt.mdp -c em.gro -o nvt.tpr
gmx mdrun -deffnm nvt -v

为了准确的得到扩散系数,你需要对数据进行充分的采样,并且最好进行多组模拟。我这里为了简单说明,因此只做了一次。

计算均方位移

根据简单扩散系数的定义,这里使用Einstein - Smoluchowski关系式来计算扩散系数:
$$D=\frac{1}{6N}\lim\limits_{t \to \infty}\frac{d}{dt}<\sum_{t}^N[r_{i}(t)-r_{i}(0)]^2>$$
计算扩散系数的步骤为先计算均方位移(MSD),然后拟合线性部分的斜率,最后除以6即为扩散系数。在GROMACS软件上使用gmx msd命令就可以得到水分子的MSD,如下: 并且最后还会输出拟合区域内的扩散系数。

1
echo -e "Water" |gmx msd -f nvt.trr -s nvt.tpr -o msd.xvg

注意:如果你不能确定拟合时间范围的话, 给出的数值最好不要直接使用。

数据分析

确定拟合区域

那么我们该如何确定计算扩散系数的拟合关联时间区域呢?方法很简单,使用MSD/6t对t作图,数据成正比的线性部分应该是一条上下稍有波动的水平线, 很容易看出来. 我们可以简单地将这段时间内扩散系数的平均值作为扩散系数。本例中的MSDMSD/6t对关联时间t作图如下:

不难看出,成线性部分的关联时间范围在1000-1800ps。

拟合计算扩散系数

确定了关联时间范围,这样得到的SPC/E水模型的扩散系数为:$2.4725*10^{-5}cm^2/s$。

实际值$2.49*10^{-5}cm^2/s$接近。

补充

对于GROMACS而言,可以通过两种方式计算自扩散系数:一种是MSD(Einstein),另一种是VACF(Green–Kubo)。前者通过gmx msd得到,后者通过gmx velaccgmx analyze得到。
有些自扩散系数是通过对速度自相关函数(VACF)来得到的,见参考资料3。这里我以0.15 mol/L NaCl溶液来计算一下钠离子和氯离子的自扩散系数。首先获得离子的VACF曲线,这里我模拟了200 ps,步长1 fs,每1 fs存储轨迹一次。最后得到了如下曲线(gmx velacc):

然后就是进行曲线积分了,采用参考资料3中的方程进行积分:
$$D_\pm=\frac{1}{3} \int_{0}^{\infty}VACF_{\pm}dt$$
支持积分的软件非常多,这里我就以xmgrace为例:

  • 绘制好上面的VACF曲线以后,依次点击Data -> Transformations -> Integration调出积分对话框,然后选择被积曲线,选择Cumulative sum,点击Accept就可以得到Sum,并且累积曲线也出来了(基线Y=0),然后Sum/3即为扩散系数,要注意自己的单位
  • 单位转换以后,我这里得到的298 k0.15 mol/L氯化钠溶液中:
    • $Na^+$的扩散系数为$1.174*10^{-5}cm^2/s$
    • $Cl^-$的扩散系数为$1.768*10^{-5}cm^2/s$
      与参考资料4吻合。

参考资料

[1] 扩散模式的分类以及扩散系数的计算
[2] A C code for calculating diffusion coefficient in MD
[3] Short communication: Investigation on dynamics and self-diffusion coefficient of [BMIM][PF6] via molecular dynamics simulations
[4] Computer simulation studies of aqueous sodium chloride solutions at 298 K and 683 K