使用Multiwfn图形化展示分子动力学模拟体系中的孔洞、自由区域

使用Multiwfn图形化展示分子动力学模拟体系中的孔洞、自由区域

文/Sobereva@北京科音  2020-Mar-5


1 前言

最近有人问我怎么获得类似下面文献中的图

此图像展示出分子动力学模拟体系中的孔洞区域,或者叫自由区域(free region)也可以,这体现出了盒子里哪些区域没有被原子占据。

之前笔者写过《使用Multiwfn可视化分子孔洞并计算孔洞体积》(http://sobereva.com/408)一文,那篇文章的做法是用于展现分子体系的特定的内部孔洞的,没法像上图这样对大量体系展示盒子里所有自由区域,而且也根本算不动。为了能得到上面的图,笔者在Multiwfn中加入了相应功能,作为主功能300的子功能1出现,在此文介绍一下。

Multiwfn可以在http://sobereva.com/multiwfn免费下载,只有2020-Mar-5及以后更新的版本才有此文的功能。不了解Multiwfn者可参看《Multiwfn FAQ》(http://sobereva.com/452)和《Multiwfn入门tips》(http://sobereva.com/167)。本文用的VMD是1.9.3版,可在http://www.ks.uiuc.edu/Research/vmd/免费下载。


2 实例

我们先直接看一个具体计算实例,之后再说相关细节。本例用到的文件是Multiwfn文件包中的examples目录下的coal.pdb,这是位网友用Lammps程序跑动力学得到的多孔介质体系的一帧,如下所示。明显能看出有不少稀疏、空白区域。

如果用文本编辑器打开此文件,会看到
CRYST1   31.064   31.100   31.093  90.00  90.00  90.00 P 1         1
其中31.064、31.100、31.093是当前这一帧的晶胞参数,即盒子的X、Y、Z尺寸,单位是埃。

启动Multiwfn,然后输入
examples\coal.pdb
300  //其它功能(Part 3)
1  //观看自由区域和计算自由区域体积
1  //设置格点并开始计算
[按回车]  //用默认的格点间距(0.3埃)。格点间距越小,要算的格点数就越多,耗时就越高,图像也会越平滑精细,因此应根据实际情况决定。0.3埃时候的图像已相当精细。

此时在屏幕上可看到以下提示
Box lengths are directly loaded from "CRYST1" field from input file
Box length in X, Y, Z:      31.064    31.100    31.093 Angstrom
Origin in X,Y,Z is           0.000     0.000     0.000 Angstrom
Grid spacing in X,Y,Z is     0.302     0.302     0.302 Angstrom
The number of points in X,Y,Z is  104  104  104   Total:     1124864
这说明程序直接从输入的pdb文件的CRYST1字段中直接读取了盒子边长。如果你用的是其它Multiwfn支持的能提供结构信息的文件,诸如xyz、mol2、gro之类的(详见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》http://sobereva.com/379),或者用的pdb文件里没有CRYST1字段,程序就会提示你手动输入盒子边长。如上还可以看到关于格点设置的一些信息。

普通四核机子花不到一分钟就能算完。算完后可以看到
Volume of entire box:   30038.649 Angstrom^3
Free volume:   14787.821 Angstrom^3, corresponding to   49.23 % of whole space
这说明整个盒子总体积,即三个边长的乘积是30038.6 Angstrom^3,自由体积是14787.8 Angstrom^3,占总体积的49.23%。自由体积是这样算的:对于每个格子,如果它距离任意一个原子的距离在此原子的范德华半径以内,这个格子就被认为是占据的。所有没有占据的格子数乘上格子体积就是自由体积。

之后看到后处理菜单,选择3,就出现了图形界面,将窗口右侧的Show molecule和Show data range都选上,可看到下图

上图中的等值面直观地展现出了盒子中未被原子占据的区域。为了得到更好的图像效果,我们选择4将格点数据导出为当前目录下的free_smooth.cub文件。启动VMD后将之载入,进入Graphics - Representation,点Create Rep,把设置改为下面的样子来显示出等值面。

最后在控制台输入pbc box和color Display Background white显示出盒子边框并改成白背景,此时图像如下

可见效果非常理想。

如果你把Isovalue设小,等值面会占更多的空间区域,如果把Isovalue设大,等值面则会收缩,且较小的等值面会消失。因此可以利用这个设置控制等值面显示的视觉效果。例如将Isovalue设为0.8后看到的图像如下,相对于上图,此时只有最显著的孔洞还可以看到,而相对次要、较小的孔洞就不显示出来了

值得顺带一提的是,以上述方式得到的孔洞的等值面与VMD的Quicksurf方式显示的分子表面是互补的。将分子的Drawing method改为QuickSurf后,VMD显示的图像如下所示

可见黄色的对应孔洞的等值面与对应分子的QuickSurf表面紧紧贴合,二者合在一起正好充满整个盒子。故这两种显示方式有高度互补性,应根据实际被研究的体系和问题恰当使用。


3 相关细节说明

下面把Multiwfn产生描述孔洞的格点数据的一些细节说一下。

在当前功能的界面里,可以看到Toggle considering periodic boundary condition选项,默认状态是Yes,即考虑周期边界条件。如果选择此选项将之状态切换为No,则计算耗时可以显著降低(约一个数量级),但构造格点数据时就不考虑原子的周期镜像了,这时候会导致误把盒子边缘的一些本来是被占据的地方显示成没有占据。下图左边是考虑周期性的,右边是没考虑周期性的,可见考虑周期性是非常重要的。

当前功能里有个选项Toggle making isosuface closed at boundary,默认状态是Yes。如果切换为No的话,等值面在盒子边缘就不是封闭的,而是镂空的,如下所示,部分镂空地方我用箭头高亮了。通常来说这样不如封闭的美观。这个设置不影响耗时。

当前功能有个选项Toggle calculating smoothed grid data of free regions,默认状态是Yes,也就是既产生“原始(primitive)格点数据”,也产生“平滑后(smoothed)的格点数据”,将之设为No不产生后者的话可以节约几分之一的时间。本文前面展示的等值面都是基于平滑后的格点数据绘制的。这里说一下两类格点数据的产生机制。原始格点数据在构造时是循环每个格点,某个点与任意一个原子的距离小于此原子的范德华半径的话,这个格点的数值就为0,否则为1,因此数值为1的格点就对应于自由区域,前面提到的自由体积也是相当于是所有数值为1的格点数的总和乘上格子体积。这种方式构造的格点数据并不适合展现孔洞区域,比如如果你在后处理菜单选择1来观看它的等值面的话,看到的是下图

可见效果巨烂,就像是奶酪中在有原子的地方挖半径等于范德华半径的窟窿。当前图像中甚至连孔洞的基本模样都看不出来。相比之下,平滑后的格点数据在构造的时候优雅得多,即每个原子被当成一个高斯函数,高斯函数半高全宽(FWHM)被设为原子的范德华半径的两倍。起初格点数据处处为1,令它减去各个原子的高斯函数分布,最后再把为负的地方设为0,就得到了Multiwfn最终给出的平滑后的格点数据了。由前面的例子可见,基于这种格点数据显示出的等值面非常光滑,而且能很好地勾勒出孔洞的位置,而且还能通过调节等值面数值来决定是只观看重要孔洞,还是也观看次要孔洞。

当前功能有个选项Set scale factor of vdW radii for calculating free volume and primitive free region,默认的scale factor是1.0。这个设的是构造前述的原始格点数据时给原子乘的范德华半径,因此影响基于原始格点数据看的等值面的效果,而且由于自由体积是对原始格点数据统计来得到的,因此这个半径也直接决定了计算出的自由体积。比如如果将前例的刻度因子设为1.5的话,算出来的自由体积就只有16.6%了。这个选项对平滑后的格点数据无影响,因为在它产生的过程中不考虑这个刻度因子。

Multiwfn的这个功能只能用于矩形盒子,如果模拟体系用的是比如三斜盒子,则这个功能不能用。此功能不仅可以用于动力学模拟得到的盒子,用于实验测定的分子晶体的晶胞等情况也完全可以。

此功能效率很高,用于甚至几万原子都可以。代码做了充分的并行化,如果体系很大算起来很慢,有以下几种方法降低耗时:
(1)用核数很多的机子算。别忘了需要将settings.ini里的nthreads设为实际的CPU的物理核数
(2)用更大的格点间距,从而减少要计算的点数。注意对于很大体系,格点间距太小的话不仅算不动,还会由于格点数太多导致内存存不下,致使程序崩溃
(3)不考虑周期性,前面已经说过了。这样做的话等值面的可靠性会下降,而且算出来的自由体积通常明显不准确