1、谈谈分子体积的计算 文 /Sobereva 2011-Sep-20 我曾在网上数次看到有人问怎么算分子体积,以及关于Gaussian 的 volume 关键词使用问题,从提问以及很多人的回答上看有不少人对这问题存在错误的认识。本文就来谈谈这个问题。 首先要认识到,分子体积不是一个可观测量,在计算方法上也不可能有唯一的定义,因此 “ 怎么计算分子体积 ” 这个问题本身就是不严格的。 分子体积就是分子表面内部空间的体积,由于分子表面也没有唯一的定义,所以不同的分子表面的定义就会给出不同的分子体积定义。首先看一个简单分子的例子 红色区域是每个原子的范德华球(以原子核为中心,半径为范德华半径的球体)的
2、叠加,这片区域就是分子的范德华体积,其表面也就被称作范德华表面。图中蓝球代表作为探针的溶剂分子(显然溶剂实际形状并不是球形,所以这个蓝球半径是 “ 等效 ” 的,在计算程序中通常是可调参数),让这个蓝球紧贴着分子范德华 表面在各处滚一遍,就产生了诸如图中绿色的轨迹,对应的表面叫做 Connolly 表面,由于溶剂分子不能触及到这个表面内的空间,所以也被叫做solvent-excluded 表面,其 内部 区 域的 体积 就 叫做solvent-excluded 体积。图中黄色是蓝球滚动时蓝球中心经过的表面,这个表面叫溶剂可及表面(其表面积就是所谓的SASA)。 最常用的体积就是范德华体积。实际
3、上范德华体积也有很多不同定义,上面介绍的原子范德华球叠加是比较简单的定义方式,通过解析几何的方法就能算出来,也可以用后面谈到的蒙特卡罗方法。这种定义存在 两个缺点: (1)原子范德华半径没有唯一定义,不同研究者给出的范德华半径存在分歧。而且对于一些金属体系没有能用的范德华半径数据。(2)没有考虑到电子效应,成键导致的电子转移、极化效果都被忽略了。 Bader 提出过一种比较严格的范德华体积的定义,消除了前述定义的弊端,并已被广泛接受,也就是若分子处于气相,则将电子密度为 0.001的等值面作为范德华表面,这种表面通常能够囊括分子 98%以上的电子密度,这种范德华体积通常比范德华球叠加得到的范德
4、华体积要大。而对于处于凝聚态的分子,考虑到分子间挤压,建议用电子密度为 0.002 的等值面作为范德华表面(显然,其体积小于0.001等值面对应的体积)。 计算范德华体积最常用的办法是蒙特卡罗方法。首先,设立一个矩形盒子,将整个分子扩住,并且各个方向上都预留一定空间以避免将范德华表面截断,记这个盒子的体积为V_box。然后,在盒子里随机分布 m 个点,依次检验这些点是否符合条件。对于范德华球叠加方式的定义,如果当前点与任何一个原子核的距离小于相应原子的范德华半径,则认为此点符合条件;而对于 Bader的定义,若当前点的电子密度大于阈值( 0.001 或 0.002),就认为符合条件。假设最后有
5、 n 个点符合条件,那么分子的范德华体积就是n/m*V_box。 如果 m 较小,那么算出来的体积是不精确的,想要增加精度,就必须增加 m。这就像人口普查,统计的人数越多结果越能反映实际国情。当然,分子越大,就需要越多的 m,才能保证平均每单位体积内随机点的数目不变,即保证精度不变。由于每次用蒙特卡罗方法计算体积时随机分布的点的位置都是不同的,因此符合条件的点数也会不同,故算出来的范德华体积的数值每次肯定会不同。然而 m 越大,结果的波动就会越小。 曾有人问什么程序计算的范德华体积更精确。这个问题太含糊,没法回答。只能概括地 说,如果想精确计算范德华体积,就需要使用 Bader的定义,用较高的
6、随机点密度去做蒙特卡罗计算,并且用比较准确的电子密度。 在 Gaussian 中,专门有个 volume关键词用于通过蒙特卡罗方法计算 Bader定义的分子范德华体积(也用来估算 Onsager溶剂模型下要用的溶质半径),同时有两个相关的 IOp。IOp(6/46=n)用来将电子密度等值面设为 n*0.0001,默认算的是 0.001 的。 IOp(6/45)可以设定每立方波尔平均有多少个随机点用于蒙特卡罗计算,默认是 20。这个数值明显偏小,结果不够精确,也导致 很多人发现每次计算体积结果都差异很大而质疑 Gaussian计算体积的功能是不是可靠。实际上,将这个数值调大 12个数量级后,体积
7、计算精度会大有改善,结果波动也会明显小很多,当然,也会更耗时。 我看到网上有人想通过让 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 程序的下一版中(目前是 2.1.2),这两种方式定义的范德华体积都可以通过蒙特卡罗方法计算。