谈谈分子体积的计算文/Sobereva 2011-Sep-20我曾在网上数次看到有人问怎么算分子体积,以及关于 Gaussian的volume关键词使用问题,从提问以及很多人的 回答上看有不少人对这问题存在错误的认识本文就来谈谈 这个问题首先要认识到,分子体积不是一个可观测量,在计算方 法上也不可能有唯一的定义,因此“怎么计算分子体积”这 个问题本身就是不严格的分子体积就是分子表面内部空间的体积,由于分子表面 也没有唯一的定义,所以不同的分子表面的定义就会给出不 同的分子体积定义首先看一个简单分子的例子红色区域是每个原子的范德华球(以原子核为中心,半 径为范德华半径的球体)的叠加,这片区域就是分子的范德 华体积,其表面也就被称作范德华表面图中蓝球代表作为 探针的溶剂分子(显然溶剂实际形状并不是球形,所以这个 蓝球半径是“等效”的,在计算程序中通常是可调参数), 让这个蓝球紧贴着分子范德华表面在各处滚一遍,就产生了 诸如图中绿色的轨迹,对应的表面叫做 Connolly 表面,由 于溶剂分子不能触及到这个表面内的空间,所以也被叫做 solvent-excluded 表 面 , 其 内 部 区 域 的 体 积 就 叫 做 solvent-excluded 体积。
图中黄色是蓝球滚动时蓝球中心经 过的表面,这个表面叫溶剂可及表面(其表面积就是所谓的 SASA)最常用的体积就是范德华体积实际上范德华体积也有 很多不同定义,上面介绍的原子范德华球叠加是比较简单的 定义方式,通过解析几何的方法就能算出来,也可以用后面 谈到的蒙特卡罗方法这种定义存在两个缺点: (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 次体积取平均 值以提高计算结果的可靠性,这是笨方法,还得写脚本去批 量执行和处理实际上直接将默认的随机点密度提高1 0 0倍, 即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),这两种方式定 义的范德华体积都可以通过蒙特卡罗方法计算。