内容发布更新时间 : 2024/5/18 8:20:48星期一 下面是文章的全部内容请认真阅读。
[图15]
与B3LYP/6-31G*波函数下的散点图相对比，在Promolecule近似下，主要区别是对应位阻效应的spike明显靠右，这说明芳香环中间的电子密度在形成分子后显著降低以减小互斥效应。并且Promolecule的散点图整体分布比较靠下。所以在观看Promolecule近似下的等值面就必须用与之前不同的等值面数值和电子密度屏蔽范围，如果仍只保留
sign(λ_2)ρ在[-0.05,0.05]区域，则最右边的spike将会被截掉很大部分；如果观看的等值面仍为0.5，则y=0.5的直线在散点图中将与一大片范围相交。所以建议保留sign(λ_2)ρ在[-0.1,0.1]区域的点，并观看RDG=0.25的等值面，这样只有这四个spike对应的区域被保留了下来，并且y=0.25的直线只与它们相交。这样得到的填色等值面如图15下方所示，色彩刻度范围设为了[-0.04,0.03]，可见和B3LYP/6-31G*下的图没有太大差别，只是芳香环中间的梭形区域变成椭圆。
虽然靠实际电子密度分析时需要保留的sign(λ_2)ρ的范围有比较通用的值[-0.05,0.05]，但是在Promolecule近似下情况比较复杂，spike的位置变化范围较大，应保留的范围不好确定，所以最好每次都考察一下散点图确定保留范围。在Multiwfn中默认的保留范围是[-0.1,0.1]，可以通过修改settings.ini里的RDGprodens_maxrho参数来改变，当sign(λ_2)ρ小于-RDGprodens_maxrho或者大于RDGprodens_maxrho的点的RDG值将被自动设为100。若此参数设为0.0，则不自动作屏蔽。
7. 用于生物大分子体系
由于在Promolecule密度下就能得到可靠结果，计算函数时速度很快，也不再需要量子化学软件计算波函数，使得这种分析方法又快又方便，能够容易地用于蛋白质、核酸等生物大分子体系。
此例将绘制DNA双螺旋结构的某一部分的RDG填色等值面图。由于这种分析方法所用的分子结构应当尽量接近平衡构型，获得这样的结构最简单的方法就是从RCSB数据库(http://www.rcsb.org/pdb/home/home.do)里面找相应的晶体结构，然后加氢。也可以自行建立DNA结构，这样序列可以自定义，可以用比如Ambertools里的NAB工具实现，之后再做能量最小化。附件里的DNA.pdb是已建好的含10个腺嘌呤-胸腺嘧啶碱基对儿的DNA片段，已通过Amber在ff99SB力场、GB溶剂模型下优化，细节见http://ambermd.org/tutorials/basic/tutorial1/。为了视觉上比较清楚，将要研究的区域是分子结构边缘区域，如下图灰色方框所示。其中黄色圆球所示的是第84、565号原子，将它们的中点作为网格中心并往外延展一定区域就能使网格涵盖感兴趣的范围。如果不清楚延展多少比较合适，可以先设定一个预期的延展范围，然后用粗糙格点计算一下，看看延展范围是否合适，再做高质量格点的计算。
[图16]
操作步骤如下
DNA.pdb //此文件已在压缩包里。假设读入某个pdb文件时提示发现未知元素，说明原子名有问题，应进行修改。 100 //功能100 1 //子功能1
16,14 //Promolecule近似下的sign(λ_2)ρ函数与RDG函数 -10 //设定延展范围 9 //延展9 bohr
7 //通过两原子中心定义网格中心
84,565 //两原子序号为84和565。注意这个序号是在Multiwfn里的序号，如果在pdb里原子序号是连贯的，比如此例，则在Multiwfn里和pdb文件里原子序号是一致的。如果pdb文件里序号不连贯，建议直接输入网格中心的坐标（网格选项中的第6项）。
120,120,120 //由于网格范围较大，所以用较多格点。
现在开始计算，计算完毕后选3将格点文件导出并用VMD作图，如图17所示。若有兴趣也可以看一下散点图，由于网
格内涉及的弱相互作用较多，所以spike比较多。
[图17]
从RDG的填色等值面图上看，上下相邻的碱基之间存在着π-π堆叠作用。腺嘌呤和胸腺嘧啶之间形成的两条氢键很清楚地通过两个蓝色圆片显现出来，在之间还存在着土黄色圆片，对应着很弱的位阻效应。红色箭头所指的是C-H---O氢键，其强度已称不上是氢键，只是范德华作用的级别。
8. 原子密度cutoff对速度的影响
在Promolecule近似下计算每个格点RDG和sign(λ_2)ρ函数时原本会计算所有原子的电子密度的贡献，当体系较大、原子较多时，比如上面例子中的体系，耗时会比较长，除非人为地将不感兴趣的区域原子删掉，但操作麻烦，还可能带来截断误差。为了进一步加快Promolecule近似下的计算，在Multiwfn中使用了原子密度cutoff的方案，对于C、H、O、N原子（由于它们最为常见），如果其电子密度在当前要算的格点位置上的数值小于1E-5，则忽视掉它，判断标准是事先推算的距离判据。这使得计算大体系时速度显著提升，对结果不会产生能查觉得到的影响。如果出于特殊原因想关掉此设定，在settings.ini里把atomdenscut设为0即可。
使用Q6600 oc 3.0G的CPU，计算中等质量格点苯酚二聚体RDG和sign(λ_2)ρ函数的总耗时： B3LYP/6-31G*波函数 312.26s Promolecule近似 10.17s
Promolecule近似+原子密度cutoff 6.16s
可见使用Promolecule近似时比使用B3LYP/6-31G*波函数时速度提升了30倍，如果基组更大，提升幅度会更大。
计算前文DNA例子的情况： Promolecule近似 753.34s
Promolecule近似+原子密度cutoff 146.48s
由于体系越大，能够少算的原子就越多，原子密度cutoff带来的加速越明显。如果不做cutoff，在计算DNA下侧的格点时还得计算所有在DNA上侧原子的贡献，白费很多时间。
9. 用RDG等值面研究成键
笔者曾在《电子定域性的图形分析》一文中介绍了拉普拉斯值函数、ELF函数、LOL函数，这三个函数对研究成键、孤对电子、原子壳层结构很有用处，但是无法用于弱相互作用分析，RDG函数分析方法是对它们重要的补充，使可视化分析能涵盖更广的范围。实际上，RDG函数对成键分析也有一定用处，从散点图上看，成键区域可看做是ρ(r)比较大的spike，在图上做横线也会与这个spike相交，因而不对ρ(r)的范围做屏蔽时也会在成键区域出现RDG等值面。具体的分析将另文讨论，这里只以典型分子吡嗪作为例子，图18上半部分是吡嗪的RDG=0.2的等值面图，下半部分是LOL=0.6的等值面图。可见，RDG函数等值面也描绘出了化学键，但是形状基本是以键轴为对称的，并且没像LOL等值面体现出共轭π键，另外RDG等值面没有描绘出氮原子的孤对电子区域。所以，RDG等值面图虽然可以用于分析成键，但并不能体现出像ELF、LOL函数那么丰富的特征，有一定局限性。
[图18]
吡嗪分子平面的RDG函数图如下所示，由此图可以想象RDG等值面轮廓随等值面数值由小到大是如何变化的，最初等值