谈谈分子体积的计算

谈谈分子体积的计算
文/Sobereva
First release: 2011-Sep-20  Last update: 2014-Jun-27


我曾在网上数次看到有人问怎么算分子体积,以及关于Gaussian的volume关键词使用问题,从提问以及很多人的回答上看有不少人对这问题存在错误的认识。本文就来谈谈这个问题。


1 分子体积的定义


首先要认识到,分子体积不是一个可观测量,在计算方法上也不可能有唯一的定义,因此“怎么计算分子体积”这个问题本身就是不严格的。

分子体积就是分子表面内部空间的体积,由于分子表面也没有唯一的定义,所以不同的分子表面的定义就会给出不同的分子体积定义。首先看一个简单分子的例子



红色区域是每个原子的范德华球(以原子核为中心,半径为范德华半径的球体)的叠加,这片区域就是分子的范德华体积,其表面也就被称作范德华表面。图中蓝球代表作为探针的溶剂分子(显然溶剂实际形状并不是球形,所以这个蓝球半径是“等效”的,在计算程序中通常是可调参数),让这个蓝球紧贴着分子范德华表面在各处滚一遍,就产生了诸如图中绿色的轨迹,对应的表面叫做Connolly表面,由于溶剂分子不能触及到这个表面内的空间,所以也被叫做solvent-excluded表面,其内部区域的体积就叫做solvent-excluded体积。图中黄色是蓝球滚动时蓝球中心经过的表面,这个表面叫溶剂可及表面(其表面积就是所谓的SASA)。


2 蒙特卡罗方法计算体积


最常用的体积就是范德华体积。实际上范德华体积也有很多不同定义,上面介绍的原子范德华球叠加是比较简单的定义方式,通过解析几何的方法就能算出来,也可以用后面谈到的蒙特卡罗方法。这种定义存在两个缺点:(1)原子范德华半径没有唯一定义,不同研究者给出的范德华半径存在分歧。而且对于一些金属体系没有能用的范德华半径数据。(2)没有考虑到电子效应,成键导致的电子转移、极化效果都被忽略了。Bader提出过一种比较严格的范德华体积的定义,消除了前述定义的弊端,并已被广泛接受,也就是若分子处于气相,则将电子密度为0.001的等值面作为范德华表面,这种表面通常能够囊括分子98%以上的电子密度,这种范德华体积通常比范德华球叠加得到的范德华体积要大。而对于处于凝聚态的分子,考虑到分子间挤压、各种形式的相互作用,会使得范德华表面有一定穿透,建议用电子密度为0.002的等值面作为范德华表面(显然其体积小于0.001等值面对应的体积)。

计算范德华体积最常用的办法之一是蒙特卡罗方法。首先,设立一个矩形盒子,将整个分子扩住,并且各个方向上都预留一定空间以避免将范德华表面截断,记这个盒子的体积为V_box。然后,在盒子里随机分布m个点,依次检验这些点是否符合条件。对于范德华球叠加方式的定义,如果当前点与任何一个原子核的距离小于相应原子的范德华半径,则认为此点符合条件;而对于Bader的定义,若当前点的电子密度大于阈值(0.001或0.002),就认为符合条件。假设最后有n个点符合条件,那么分子的范德华体积就是n/m*V_box。

如果m较小,那么算出来的体积是不精确的,想要增加精度,就必须增加m。这就像人口普查,统计的人数越多结果越能反映实际国情。当然,分子越大,就需要越多的m,才能保证平均每单位体积内随机点的数目不变,即保证精度不变。由于每次用蒙特卡罗方法计算体积时随机分布的点的位置都是不同的,因此符合条件的点数也会不同,故算出来的范德华体积的数值每次肯定会不同。然而m越大,结果的波动就会越小。

曾有人问什么程序计算的范德华体积更精确。这个问题太含糊,没法回答。只能概括地说,如果想精确计算范德华体积,就需要使用Bader的定义,用较高的随机点密度去做蒙特卡罗计算,并且用比较准确的电子密度。

在Gaussian中,专门有个volume关键词用于通过蒙特卡罗方法计算Bader定义的分子范德华体积(也用来估算Onsager溶剂模型下要用的溶质半径),同时有两个相关的IOp。IOp(6/46=n)用来将电子密度等值面设为n*0.0001,默认算的是0.001的。IOp(6/45)可以设定每立方波尔平均有多少个随机点用于蒙特卡罗计算,默认是20。这个数值明显偏小,结果不够精确,也导致很多人发现每次计算体积结果都差异很大而质疑Gaussian计算体积的功能是不是可靠。实际上,将这个数值调大1~2个数量级后,体积计算精度会大有改善,结果波动也会明显小很多,当然,也会更耗时。

我看到网上有人想通过让Gaussian计算100次体积取平均值以提高计算结果的可靠性,这是非常笨的行为,还得写脚本去批量执行和处理。实际上直接将默认的随机点密度提高100倍,即IOp(6/45=2000),算一次就行了,完全是等价的。

用Volume关键词时有两点需要注意。第一点是对于后HF方法不要忘了加density关键词,否则密度是HF下的,但实际上对结果的影响很小,一般还不如蒙特卡罗方法的随机性的影响大。另外就是输出的结果单位问题,会有诸如这样的输出
 Molar volume =  532.998 bohr**3/mol ( 47.564 cm**3/mol)
这里bohr**3/mol应为bohr**3。另外括号里的数据只是通过简单单位换算得到的,而并非是通过什么经验公式得到的能和实际物质密度相比较的结果。

用Bader的定义计算范德华体积,由于涉及到计算很多点的密度,所以比使用范德华球叠加方式的定义要慢很多,对于大分子,则Bader的定义会因为太昂贵而无法使用。在笔者的Multiwfn程序中,这两种方式定义的范德华体积都可以通过蒙特卡罗方法计算。此程序可以在http://sobereva.com/multiwfn免费下载。

在Multiwfn程序的settings.ini文件里,如果MCvolmethod被设为了1,代表使用范德华球叠加方式定义体积;如果被设为了2(默认),代表使用Bader的定义。假设要计算水分子的范德华体积,那么启动程序后首先载入水分子的.wfn文件或者.fch文件,然后选100,再选3。程序会提示你输入i、x和k值(如果是用范德华球叠加方式,只需要输入i)。这里i代表将会有100 * 2^i个随机点分布在盒子内,x代表等值面的电子密度数值,k代表盒子在分子周围预留的空间为k乘以预置的原子范德华半径以避免等值面被盒子截断。通常情况下,输入9,0.001,1.7就可以了,对于较大的分子,建议将i值设大来得到更准确的结果。输入i、x和k后程序会立刻给出结果:
The molecular volume:  171.982 bohr^3, (   25.485 angstrom^3,   15.347 cm^3/mol)
这个结果与Gaussian的Volume关键词的结果可以相互对照。


3 基于MT算法计算体积


另外,Multiwfn还有一种计算体积的方式,就是通过主功能12(定量分子表面分析)来计算,在输出分子表面信息的同时也会一起把体积信息输出出来,这也是Bader定义的范德华体积,但是用的不是蒙特卡罗算法计算而是用本人提出的改进的Marching Tetrahedron (MT)算法,算法细节见J. Mol. Graph. Model., 38, 314 (2012)。你会发现这两种算法得到的体积结果相符很好。原理上,只要蒙特卡罗用的点数无穷大,MT算法用的电子密度格点精度无穷高,两种方法的结果是完全相同的。我个人更倾向于用MT方法来计算范德华体积,因为没有蒙特卡罗方法的那种随机性。

MT算法简单来说就是一种构建等值面的方法,需要格点数据进行输入。很多可视化软件读取cube文件后显示的等值面实际上就是用MT方法产生的。Multiwfn在做定量分子表面分析过程中会先生成电子密度格点数据,然后由此依靠MT算法构建电子密度为0.001等值面(也可以是0.002等,用户随意设定),分析这个面上的各种性质,同时输出这个面里面包含的空间大小,因此就是Bader定义的范德华体积了。

关于Multiwfn的定量分子表面分析功能详见《使用Multiwfn的定量分子表面分析功能预测反应位点、分析分子间相互作用》(http://sobereva.com/159)、手册3.15节的原理介绍以及手册的4.12节的例子。这里仅是举一个最简单的例子说明怎么用这个功能计算范德华体积。启动Multiwfn后,依次输入
examples\phenol.wfn        //体系的波函数文件,可由Gaussian等程序得到,具体方法见Multiwfn手册第四章开头
12  //主功能12
2   //切换映射的函数
-1  //用户自定义函数
0  //开始运算
算完后就会看到输出的范德华体积
Volume:   833.28906 Bohr^3  ( 123.48073 Angstrom^3)
还有其它一些信息,包括范德华表面积
Overall surface area:         472.15159 Bohr^2  ( 132.21593 Angstrom^2)
之所以在选0之前先把映射函数切换为了“用户自定义函数”,是因为默认会对分子表面的静电势进行分析,但是静电势计算很耗时,而我们当前只需要得到体积而不需要分析静电势,所以把映射的函数设为了在默认情况下不需要任何耗时的“用户自定义函数”。默认情况下分子表面的定义用的是电子密度0.001等值面,比如想改成0.002等值面来得到凝聚相下分子范德华体积,在选0之前先选1设定表面定义,选择电子密度等值面,然后输入0.002即可。


4 其它问题


有人在某些数据库上查到“分子体积”或“范德华体积”,这里要强调的是,不要盲目地相信这些查到的数据!很多这些数据库里的数据仅仅是通过非常简单的模型(比如不同元素的原子数和连接关系)估算出来的,显然不如通过电子密度等值面算出来的体积更有合理性。使用、讨论这些查到的数据之前,一定先搞清楚这数据是怎么得来的!

虽然范德华体积没有唯一的定义,但在实验上一种相对有意义的范德华体积的定义是从物质的密度和分子质量上反推出来的分子体积(下面简称“实验范德华体积”。注意通过其它实验手段也可以定义不同的实验范德华体积)。有人会发觉并且抱有疑问,怎么这种实验范德华体积和0.002等值面包围的体积数值差距很大?其原因也是很明显的。前文计算分子体积得到的是对单个分子的描述,是单分子的性质。而密度是物质的宏观性质,由此反推出来的实验范德华体积本质上也是宏观性质。分子间的各种各样复杂的相互作用直接影响物质的宏观性质(范德华作用、静电作用、分子间挤压、构象变化等),也因此影响实验范德华体积的因素十分复杂,原理上就不可能仅仅基于单个分子电子结构性质就准确得到这种实验范德华体积,除非去做分子动力学模拟来获得密度并反推出体积。但是,通过前面的方法计算出的分子范德华体积,必然和实验范德华体积有相关性,并且对于同一类物质,偏差会是较为系统的。因此,可以自行拟合一个方程(通常线性的就可以),来将二者关联起来。

添加新评论