使用Multiwfn计算分子的动力学直径

sobereva

使用Multiwfn计算分子的动力学直径

Using Multiwfn to calculate dynamic diameter of molecules

文/Sobereva@北京科音

First release: 2019-Aug-19  Last update: 2024-May-4

1 相关知识&计算原理

分子直径有很多不同定义。在分子分离研究领域中经常用动力学直径（kinetic diameter）体现分子尺寸，这个量没有一般性的方法，直到现在，主要来源还是Breck很老的1974年的Zeolite Molecular Sieves; Structure, Chemistry and Use一书，数值是通过气体在不同温度和压力下的第二维里系数实验和Lennard-Jones势的形式推出来的，动力学直径对应LJ作用为0的情况。

在J. Phys. Chem. A, 118, 1150 (2014)文中，作者提出了基于分子电子密度等值面的估计分子动力学直径的方法。对于文中测试的一系列小分子，计算结果和被普遍采用的Breck书里的那些值有不错的线性相关性，因此此文的做法可以作为一个普适型的算动力学直径的做法。

分子表面可以通过电子密度等值面来定义，这点我在很多博文里提过，比如《使用Multiwfn和VMD计算分子表面积和片段表面积》（http://sobereva.com/487）、《谈谈分子体积的计算》（http://sobereva.com/102）、《使用Multiwfn的定量分子表面分析功能预测反应位点、分析分子间相互作用》（http://sobereva.com/159）等。JPCA文中的思想是以这种方式定义分子表面，然后测量最短方向的等值面两端的距离，即下面这些例子中的箭头夹着的距离。

前述JPCA文中在PBE0/def2-TZVP级别下考察了不同电子密度下算的动力学直径与Breck的值的线性关系，发现用0.0015 a.u.作为等值面定义效果最好。如下所示，这样计算的几个小分子的动力学直径（纵坐标）与Breck书中的动力学直径（横坐标）的线性关系不错。

如JPCA文中表2所示，用0.0015 a.u.作为等值面数值计算后，原理上最好再除以拟合出来的系数1.025。

此文的方法思想非常简单，也没有考虑分子在实际环境中电子的可极化效应导致对动力学直径产生的影响。但不管怎么说，这个方法还是比较有实用性的，如果被考察的分子查不到现成的动力学直径，不妨用这个方法粗略估计一下。

JPCA这篇文章的方法可以通过Multiwfn的定量分子表面分析功能非常容易地实现，下面就以CO2分子为例讲一下如何计算。对Multiwfn完全不了解者参看《Multiwfn FAQ》（http://sobereva.com/452）、《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）。使用Multiwfn以这种方式计算动力学直径已经被很多文章所使用，例如Nature (2024) DOI: 10.1038/s41586-024-07342-y。


2 计算实例：CO2

这里我们首先用Gaussian程序优化CO2并产生它的波函数。笔者用的是Gaussian16 A.03，计算级别用的和JPCA里相同，即PBE0/def2-TZVP。输入文件如下

%chk=C:\CO2.chk
# PBE1PBE/def2TZVP opt
[空行]
Title Card Required
[空行]
0 1
C                  0.00000000    0.00000000    0.00000000
O                  0.00000000    0.00000000    1.25840000
O                  0.00000000    0.00000000   -1.25840000

将CO2.chk用formchk转换成CO2.fch，之后就可以用Multiwfn计算了，此文件可以从这里下载：http://sobereva.com/attach/503/file.rar。用Multiwfn载入CO2.fch，然后输入以下命令
12  //定量分子表面分析
1   //选择表面的定义
1   //用电子密度等值面作为表面
0.0015  //用0.0015 a.u.作为等值面数值
6  //开始分析（但不考虑映射的函数，因为这不是我们目前感兴趣的）
6  //将表面顶点坐标导出为当前目录下的vtx.pdb

现在，通过vtx.pdb我们就可以考察动力学直径了，有两种方法可以实现。下文的“方法1”最普适，可以用于任意形状体系，但需要借助VMD，且眼力不能太差，要有耐心；而“方法2”的做法省事得多，不过没那么灵活，不适合形状复杂的情况。


方法1：通过VMD测量表面顶点间距离

在http://www.ks.uiuc.edu/Research/vmd/下载VMD并安装，启动后将vtx.pdb拖入其中，在Graphics - Representation里将Drawing Method设为Points。然后在VMD Main窗口选择Display - Orthographic使用正交视角以便于选取顶点。将图形窗口拉大，按键盘上的2进入距离测量模式，恰当旋转分子使得要选的顶点处在容易被选中的位置，然后按照JPCA文中的示意图，点击能反映动力学直径的两个位置的表面顶点。恰当选择后，该体系两个角度的图像如下所示

从VMD的文本窗口中可以读到精确数值：
Info) Added new Bonds label MOL1:C/MOL1:C = 3.385527
即动力学直径是3.385埃，和Breck书中的3.31埃相符很好，和JPCA文中表1中的3.469埃也基本吻合（JPCA这篇文章没说具体怎么测量的，但至少我们的测量方法是肯定恰当的）。

用此方法选取顶点时往往会碰到一个问题，就是从某个角度看选取的顶点似乎很合适，但换一个角度看又不那么合适。选取的技巧只可意会不可言传，试多了自然就懂了。点上去之后如果发现不合适想删除标签，可以进入Graphics - Labels，在Atoms和Bonds里圈上条目后点Delete即可删除。


方法2：通过Multiwfn直接统计表面顶点的最大、最小坐标

用这个方法前，必须先保证当前体系的最短方向冲着X、Y、Z坐标轴之一。比如对于CO2.fch，利用Multiwfn主功能0一看就可以看出最短轴是X或Y方向，如下所示，因此这个体系可以用下面的方法考察动力学直径。

启动Multiwfn并载入vtx.pdb，进入主功能100，然后进入子功能21，这个功能专门考察体系的几何信息。然后输入all，从输出文件中可以找到以下信息
Minimum X is   -1.69400000 Angstrom, at atom   134(C )
Minimum Y is   -1.69400000 Angstrom, at atom  2112(C )
Minimum Z is   -2.75500000 Angstrom, at atom  2177(C )
Maximum X is    1.69300000 Angstrom, at atom  4249(C )
Maximum Y is    1.69400000 Angstrom, at atom  2244(C )
Maximum Z is    2.75600000 Angstrom, at atom  2179(C )

由于Y轴方向原子位置（当前语境下对应表面顶点位置）最小值是-1.694，最大值是1.694，因此CO2的动力学直径是1.694*2=3.388埃。此数值和“方法1”肉眼测量的3.385埃很接近，显然用这种方法比方法1省事得多，还可以避免肉眼测量可能造成的不准确。但实际研究的分子往往比此例复杂得多，此方法未必总适用。有时虽然可以用，但Gaussian算之前可能需要将分子先恰当旋转到合适朝向，并且计算时加上nosymm关键词避免自动被旋转，详见《谈谈Gaussian中的对称性与nosymm关键词的使用》（http://sobereva.com/297）。

值得一提的是，利用Multiwfn和VMD，不仅可以像上文这样考察分子的外径，还能通过测量表面顶点距离的方式考察环状、笼状分子的内径。在《一篇最全面、系统的研究新颖独特的18碳环的理论文章》（http://sobereva.com/524）介绍的笔者的论文中就用Multiwfn+VMD测量了18碳环的内径，并讨论什么样的分子有可能穿过去，感兴趣的读者建议看看。

Multiwfn还专门有计算分子孔洞直径的功能，见《使用Multiwfn计算分子和晶体中孔洞的直径》（http://sobereva.com/643），只不过这是基于原子坐标和范德华半径，而非基于电子密度等值面算的。

youyno

大神棒棒哒！

Shine剪水

http://sobereva.com/190好像和这篇博文的方法3,4很类似

ydba

sob老师，可以通过http://sobereva.com/426这篇帖子中提到的“将分子摆正”后的new.pdb文件，再读取表面顶点的最大、最小坐标吗得到动力学直径吗？

sobereva

ydba 发表于 2023-7-7 05:07
sob老师，可以通过http://sobereva.com/426这篇帖子中提到的“将分子摆正”后的new.pdb文件，再读取表面顶 ...

你可以用Multiwfn摆正后的结构计算波函数文件，按本文做法计算表面顶点位置

zhangweiwei

卢老师好！@sobereva 您说：“如JPCA文中表2所示，用0.0015 a.u.作为等值面数值计算后，原理上最好再乘上拟合出来的系数1.025。”，文中得到的拟合曲线斜率是1.025，如图。y是QM结果，x是breck实际结果。应该是除以1.025才对吧？

sobereva

zhangweiwei 发表于 2024-5-4 00:22
卢老师好！@sobereva 您说：“如JPCA文中表2所示，用0.0015 a.u.作为等值面数值计算后，原理上最好再乘上拟 ...

除以

zhangweiwei

卢老师您好！@sobereva 这个方法我已经基本掌握，并验算了CO2。但是我按照您这个帖子，尝试验证了一下CH4的计算，发现有较大偏差（算出来4.537，原文4.046）。请问是我步骤不对，还是没有领悟原文/帖子的精髓？在第二步一定需要fch文件吗？因为如果我旋转了，就只能输出pdb/xyz文件，然后用pdb/xyz文件计算0.0015的等势面。
我的步骤如下：
gausssian view新建一个CH4.gjf，如帖子中一样操作，计算出chk转成fch。由于此时C-H键并不是平行于xyz轴，于是按照【Multiwfn中非常实用的几何操作和坐标变换功能介绍http://bbs.keinsci.com/forum.php ... 24674&fromuid=62697】，将C-H键旋转与z轴平行，然后生成pdb/xyz文件。然后将这个pdb/xyz文件按照您帖子里的步骤计算生成vtx.pdb得到如下信息：
Minimum X is   -2.20800000 Angstrom, at atom      40(C )
Minimum Y is   -2.20800000 Angstrom, at atom    1710(C )
Minimum Z is   -1.99600000 Angstrom, at atom    2491(C )
Maximum X is    2.44900000 Angstrom, at atom    4995(C )
Maximum Y is    2.44900000 Angstrom, at atom    2302(C )
Maximum Z is    2.54100000 Angstrom, at atom    2493(C )
可以看到最短方向Z轴，长度为4.537，与原文4.046相差很大。
文件附在后面。如果是您，您会怎么做？恳请解答，谢谢！

sobereva

zhangweiwei 发表于 2024-5-5 05:14
卢老师您好！@sobereva 这个方法我已经基本掌握，并验算了CO2。但是我按照您这个帖子，尝试验证了一下CH4的 ...

Multiwfn里旋转后，导出gjf文件，Gaussian重新算个单点（加上nosymm避免Gaussian自动修改朝向）得到新的波函数文件，重新再按照博文里的步骤计算直径就完了

sobereva

今天看到一篇Nature文章用Multiwfn通过本文介绍的方法计算了动力学直径
Xiangjie Zhang, Tao Yan, Huaming Hou, et al., Regioselective hydroformylation of propene catalysed by rhodium-zeolite, Nature (2024) https://www.nature.com/articles/s41586-024-07342-y

chen0201

sobereva 发表于 2024-5-6 12:51
今天看到一篇Nature文章用Multiwfn通过本文介绍的方法计算了动力学直径
Xiangjie Zhang, Tao Yan, Huaming ...

请问老师，我按照您的方法计算了分子动力学直径，但是对于苯乙炔这种分子来说，动力学直径是用躺着的苯环（分子动力学直径小~4）还是用立着的苯环（~6.5）来描述比较合适呢？或者有更好定量描述类似苯乙炔分子空间位阻效应的方法吗？

sobereva

chen0201 发表于 2024-7-18 15:19
请问老师，我按照您的方法计算了分子动力学直径，但是对于苯乙炔这种分子来说，动力学直径是用躺着的苯环 ...

躺着的
文中CO2的例子体现出各个方向跨度不同时怎么考虑