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