本帖最后由 Uus/pMeC6H4-/キ 于 2025-7-8 11:09 编辑
用VMD看Gaussian隐式溶剂模型(IEFPCM/SMD)计算的溶质孔洞及定性表面电荷分布
Visualizing solute cavity and qualitative surface charge distribution in Gaussian implicit solvation model (IEFPCM/SMD) calculations with VMD Uus_pMeC6H4- (\ Init post 2025-07-05 /) Last edit 2025-07-08
对溶液相物质做量子化学计算时，溶剂对溶质的极化是隐式溶剂模型考虑的重点。可极化连续介质(PCM)类隐式溶剂模型需要构建溶质孔洞、将孔洞表面划分为大量面积微元碎片、求解每块小碎片上的显著表面电荷。对于广泛应用的IEFPCM和SMD而言，溶质孔洞通常由以原子坐标为中心、原子范德华半径与缩放因子的乘积为半径的球叠加而成，但也有溶剂排斥表面SES、溶剂可及表面SAS等选择，可以参考此帖的原理说明（注意此帖是近8年前写的，软件运行排障的内容未必仍然适用）。
对Gaussian计算而言，用GaussView 6打开输出文件，可以在Results - PCM Solvation Cavity查看溶剂孔洞，调节选项包括Solid、Mesh、Point等显示风格和单一、分原子等着色方式，如这张科音基础量子化学培训班的幻灯片所示（来自此帖）。
本帖将提供一个读取Gaussian计算输出的文件、在VMD里观看溶质孔洞的方法，以获得调节灵活度更高的效果。相关脚本在下述压缩包里，下载后需将其中solcav.vmd和GeomView.vmd两个文件移到VMD目录下，修改vmd.rc（或.vmdrc）文件添加source solcav.vmd一行，保存并重新启动VMD。在VMD命令行里输入solcav --help命令，如果产生一长串帮助信息即说明脚本生效。
solcav.zip
(44.25 KB, 下载次数 Times of downloads: 11)
先以水中乙醇胺盐酸盐HOCH2CH2NH3+···Cl-以氢键联系的紧密离子对直链构象（忽略其他参与形成氢键的第一配位层水分子）为例说明下脚本的用法。作为对比，可以先用http://sobereva.com/443的方法快速画一个静电势着色的分子范德华表面及静电势极值点图（turbo渐变色方案从蓝到红对应-0.10到0.10 a.u.的静电势；原理上溶液这样的凝聚相用电子密度0.002 a.u.等值面更合适，不过这里姑且还用0.001 a.u.的等值面；图片用Tachyon渲染，参考http://sobereva.com/624给碳用2号gray、氯用19号green2的颜色，下同）。
将下述输入文件hoch2ch2nh3cl.gjf里的{empty}换成空行后，提交任务做一个单点计算，只需1 min就能完成
- #p wb97xd/def2tzvpp scrf=(solvent=water,read)
- {empty}
- hoch2ch2nh3cl
- {empty}
- 0 1
- C -2.22299271 0.35589912 0.00000191
- H -2.39627414 0.97423192 -0.88671406
- H -2.39626921 0.97422672 0.88672248
- C -0.79755136 -0.15874777 -0.00000395
- H -0.60337023 -0.76121648 0.88397493
- H -0.60337422 -0.76120755 -0.88398980
- N 0.17672744 0.95827490 -0.00000050
- H 0.06924251 1.54910878 0.82045304
- H 1.16376407 0.57450768 -0.00000812
- H 0.06923490 1.54912074 -0.82044441
- O -3.05552475 -0.78168174 0.00000093
- H -3.96913577 -0.49658528 0.00000088
- Cl 2.94097367 -0.30820971 0.00000079
- {empty}
- PrintSpheres
- GeomView
- {empty}
- {empty}
此时目录下除了主要的输出文件hoch2ch2nh3cl.out，还多产生了两个文件charge.off和points.off。points.off不用管，将charge.off重命名为hoch2ch2nh3cl.off，然后把hoch2ch2nh3cl.out和hoch2ch2nh3cl.off复制到VMD文件夹里，启动VMD并在其命令行输入solcav hoch2ch2nh3cl得到下述效果，即为溶质孔洞图形。
若在命令末尾加个1，即solcav hoch2ch2nh3cl 1这样（上面尝试完以后需先重启VMD再操作），则得到下述效果，溶质孔洞表面出现很多小圆球。
那么solcav.vmd是如何工作的？这就先从上面Gaussian的输入文件讲起。以scrf关键词启动隐式溶剂模型时，可以用scrf(read)选项在末尾添加额外输入，包括http://sobereva.com/327里提到的各种自定义选项，具体看文档。此处添加的PrintSpheres选项确保Gaussian会在主要输出文件里打印溶质孔洞的具体信息（没有PrintSpheres时仅<=20个球的体系会默认打印），包括球心坐标(Xe,Ye,Ze)、初始原子半径Re0、缩放因子Alpha、归属原子on等，供solcav.vmd读取并产生VMD的结构。Gaussian默认给IEFPCM用UFF力场定义的范德华半径和1.1的缩放因子，而给SMD用分元素定义的范德华半径（具体参考http://sobereva.com/613第一段）和1.0的缩放因子。
另一边，GeomView是一个三维可视化软件的名称，有这个选项时Gaussian会以该软件自定义的格式输出charge.off和points.off两个文本文件。观察发现：
(1)这两个文件包含溶质孔洞表面各面积微元碎片对应的顶点的坐标，与输出文件对照可知为Lebedev-Laikov方法生成的积分格点（参考http://sobereva.com/69第3.2节）、总数在GePol: Number of points一行给出，默认密度每平方埃5.0个格点，可以用scrf(read)额外输入的PDens选项调节；
(2)charge.off与points.off均给每个顶点定义了颜色的RGBA值，其中charge.off的颜色是由显著表面电荷决定的，具体来说负为蓝、零为白、正为红，越接近最小值（最大值）则蓝色（红色）越深，而蓝白红渐变色中RGB只有15个离散的值而不是连续变化的，按我前一帖用matplotlib或gnuplot给VMD定义渐变色标尺的方法用make_palette.py展示如下。
基于这两点，solcav.vmd载入.off文件时将每个顶点当作一个小氢原子，其charge属性根据颜色定为-7到7的一个整数，以便用Charge的着色方式显示。solcav.vmd也会相应自动加载GeomView.vmd，修改33到1056号颜色的RGB值，处理完毕打开Graphics - Colors窗口切换到Color Scale栏即可看到上述渐变色方案。（不知为何，在我测试的几个体系中总有不属于上述15个RGB值的颜色出现；solcav.vmd考虑这些outlier时仍然加载但隐藏显示，具体隐藏了哪些看显示方式的选区。）修改渐变色产生的兼容问题在我的定义渐变色标尺一帖末尾有详细解释，此处不重复。
对照静电势着色的分子表面图也能确认显著表面电荷的分布，比如上面的HOCH2CH2NH3+···Cl-中，氧原子和氯原子附近的静电势都比较负，都会在溶质孔洞表面极化产生正的显著表面电荷，但氯离子比羟基氧带的负电荷更多，所以氯周围的显著表面电荷更正，红色更深。然而讨论仅止于此，毕竟从charge.off各点的颜色倒推的显著表面电荷只分了离散的15个档次、只能定性体现同一体系内的相对大小，无法像multiwfn的定量分子表面分析功能（参考http://sobereva.com/196）一样给出具体的数值，也无法用于比较不同体系。本帖的标题特意用了“定性”一词就是为强调这一点。
最后，solcav.vmd除载入上面溶剂孔洞及表面顶点外，也会把输出文件里体系本身的原子笛卡尔坐标给读取了，默认先找标准朝向（或Z-matrix朝向），如果找不到就看输入朝向，而二者的区别在http://sobereva.com/297有详细讲解。原子坐标、溶剂孔洞、表面顶点是三个独立的molecule，均用原子（而不是如draw sphere画的球）来展示，方便按需调整原子选区和显示风格。强烈建议用Tachyon或POV-Ray等渲染图像，效果比图形窗口里直接看到的或者Snapshot截屏的好看多了。
| 
发表于 Post on 2025-7-5 02:18:56
|
只看该作者 Only view this author
|只看大图 View big picture
|倒序浏览 Reverse view
|阅读模式 Reading model
收藏 Add to favorites
楼主 Author