扩散系数计算

本文主要简介了通过GROMACS软件计算扩散系数的步骤。并提供一些不同方向扩散系数实例步骤。

体系构建和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作图脚本以及图如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#!/bin/bash
awk '
BEGIN {
HeadXvg()
}
/^[^#@]/ {
if($1+0 != 0) {
print $1, $2/(6*$1)
}
}

function HeadXvg() {
print "@ xaxis label \"Time (ps)\""
print "@ yaxis label \"MSD/(6*t) (nm\\S2\\N/ps)\""
print "@TYPE xy"
}
' msd.xvg > msd_6t.xvg

不难看出,成线性部分的关联时间范围在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吻合。

补充[2020.03]

最近看到有人问如何做【如何计算液体的分层扩散系数?】。引用分子模拟周刊:第 12 期中的点评:

这个问题有年头了, 问过的人也不在少数, 但到目前也没有太好的解决方法, 因为每个粒子不是固定于某层的, 所以要处理只能采用类似计算滞留时间的方法, 比较麻烦. 文献上还有其他的一些处理方法, 看起来都比较复杂. 我一直在找机会把这个问题处理一下, 可大多数问这个问题的人也就只是问问, 并不见得真是铁了心要把它解决掉. 既然如此, 我又何必越粗代庖呢?

我对这方面也是新手,不过正如上面点评所述,看了很多文献上的做法,单看给出的计算公式就够你喝一壶的了,最常用的是先求MFPT和自由能,然后计算不同位置的扩散系数【文献5】。我也是懵逼中…
不过后来又想了想,既然要分区域,又要保证分子或者原子在改区域内,那就得保证取样时间差合适。何不用有限差分的思想来给出一种使用于GROMACS模拟后期分析的思路呢?
经过一番考虑,我觉得结合GROMACS自带的一些分析工具gmx trjconvgmx selectgmx msd,然后写一点自动化代码应该是可以勉强实现一下的。下面给出我的计算思路,对不对也不太清楚,但感觉方向大致是对的吧!


给出一个模拟体系,中间为石墨烯的TIP3P水体系。即在Z轴方向被石墨烯片平分,Z轴方向体系长度为4 nm
实现代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#!/bin/bash
tpr="md.tpr"
xtc="md.trr"
gro="conf.gro"
ndx="grondx.ndx"

for z in `seq 0 0.01 3.99`;do
for(( i=0;i<=50;i=i+3 ));do
let t=i+3
echo 0 |gmx trjconv -f $xtc -s $tpr -b $i -e $i -o $gro > /dev/null 2>&1
zl=`awk -v z=$z 'BEGIN{print z+0.01}'`
gmx select -s $gro -on $ndx -select "resname SOL and z>$z and z<=$zl" > /dev/null 2>&1
gmx msd -f $xtc -s $gro -n $ndx -b $i -e $t -o msd_${i}.xvg > /dev/null 2>&1
rm *#
done
echo $z
grep "cm^2/s" msd_*.xvg| awk '{sum += $(NF-4)} END{print sum/NR}' >> mmsd.txt
done
rm msd_* conf.gro

为了节省时间,上面的代码取样并不充分,但也能够说明问题。因为GROMACS自带的gmx density能够得到沿Z轴方向的水分子密度,可以和你计算的扩散系数进行对应比对。因此最终我得到了如下图:

实际的纯TIP3P水模型扩散系数为$5.19*10^{-5}cm^2/s$,是不是看上去还挺不错的?石墨烯存在的中间位置水分子扩散系数明显要低,也基本上呈现对称性。

行吧,,,就这那样吧。纯属好玩~

补充[2020.05]

最近又看到一篇关于如何计算不同方向上的扩散系数文章【文献5】。也是采用MFPT方法计算的,但我觉得这个方法要比以前的文章简单的多,特别是简化了一些步骤,使得自己也可以通过简单的编程来进行统计计算。
比如你想计算某个界面垂直Z方向沿路径的水分子扩散系数,主要步骤如下:

  • 1、计算Z方向粒子密度分布,得到沿该方向的自由能曲线$\beta U_z$

  • 2、写程序计算水分子氧原子在Z方向区间$\pm{\delta_z\over2}$内粒子平均首次通过边界$[Z-L, Z+L]$的时间$<\tau_z>$,以及通过左右边界的可能性$P_+z$和$P_-z$

  • 3、求其$U_z$的一阶导数$U’_z$,得到变量$h$ = $\beta U’_z$;计算$\theta (h, \kappa)$; 计算Z方向扩散系数$D_z$

这里涉及到第二步需要自行提取轨迹中水分子氧原子的坐标,然后进行计算统计。这里我给出一部分代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
//略
double incr = 0.03, deltaz = 0.02, L = 0.1, z0 = 0.1, z = 0, sum; //nm
real *temp;
int *ok, count, lcount, i;
gmx_rmpbc_t gpbc;

for(int j = 0; j < 134; j++) {
z = z0 + j*incr;
natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
gpbc = NULL;
gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);

sum = 0;
count = lcount = 0;
snew(temp, isize[0]);
snew(ok, isize[0]);

do {
gmx_rmpbc(gpbc, natoms, box, x);
for(i = 0; i < isize[0]; i++) {
if(ok[i] == 0 && x[index[0][i]][2] >= z-deltaz/2 && x[index[0][i]][2] <= z+deltaz/2) {
temp[i] = t;
ok[i] = 1;
}
else if(ok[i] == 1 && (x[index[0][i]][2] <= z-L || x[index[0][i]][2] >= z+L)) {
sum += t-temp[i];
ok[i] = 0;
count++;
if(x[index[0][i]][2] <= z-L) {lcount++;}
}
}
}
while (read_next_x(oenv, status, &t, x, box));
sfree(temp);
sfree(ok);
fprintf(fmftp, "%8.3lf %8.6lf %8.3lf %8.3lf\n", z, sum/count,
double(lcount)/count, 1-double(lcount)/count);
fflush(fmftp);
}
//略

这样我们就可以得到$<\tau_z>$,$P_+z$和$P_-z$这几个参数了。最后统计一些数据可以得到类似下面的扩散系数$D_z$随Z方向距离(垂直于界面)的距离变化图了。

参考资料

[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
[5] First-passage fingerprints of water diffusion near glutamine surfaces