分子体积的计算.doc

谈谈分子体积的计算文/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次体积取平均值以提高计算结果的可靠性，这是笨方法，还得写脚本去批量执行和处理。实际上直接将默认的随机点密度提高100倍，即IOp(6/45=2000)，算一次就行了，完全是等价的。 用Volume关键词时有两点需要注意。第一点是对于后HF方法不要忘了加density关键词，否则密度是HF下的，但实际上对结果的影响很小，一般还不如蒙特卡罗方法的随机性的影响大。另外就是输出的结果单位问题，会有诸如这样的输