用VMD看Gaussian隐式溶剂模型(IEFPCM/SMD)计算的溶质孔洞及定性表面电荷分布

Uus/pMeC6H4-/キ

用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: 23)

2025-7-8 11:08 上传 Uploaded

点击下载
Click to download

先以水中乙醇胺盐酸盐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就能完成

1. #p wb97xd/def2tzvpp scrf=(solvent=water,read)
   
2. {empty}
   
3. hoch2ch2nh3cl
   
4. {empty}
   
5. 0 1
   
6. C                 -2.22299271    0.35589912    0.00000191
   
7. H                 -2.39627414    0.97423192   -0.88671406
   
8. H                 -2.39626921    0.97422672    0.88672248
   
9. C                 -0.79755136   -0.15874777   -0.00000395
   
10. H                 -0.60337023   -0.76121648    0.88397493
    
11. H                 -0.60337422   -0.76120755   -0.88398980
    
12. N                  0.17672744    0.95827490   -0.00000050
    
13. H                  0.06924251    1.54910878    0.82045304
    
14. H                  1.16376407    0.57450768   -0.00000812
    
15. H                  0.06923490    1.54912074   -0.82044441
    
16. O                 -3.05552475   -0.78168174    0.00000093
    
17. H                 -3.96913577   -0.49658528    0.00000088
    
18. Cl                 2.94097367   -0.30820971    0.00000079
    
19. {empty}
    
20. PrintSpheres
    
21. GeomView
    
22. {empty}
    
23. {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截屏的好看多了。

Uus/pMeC6H4-/キ

刚做了个补充材料，把几种不同情况的溶质孔洞及表面顶点导出为Wavefront OBJ格式并插入.pptx幻灯片中。原始文件尺寸达711 MB，但发现用xz适当压缩并用split切割以后能满足论坛“单个文件不超过9.8 MB，一日文件总大小不超过48.8 MB”的要求，故上传附件如下。
Supplementary_Information.tar.xz.part-00 (7 MB, 下载次数 Times of downloads: 15)

2025-7-6 02:12 上传 Uploaded

点击下载
Click to download
Supplementary_Information.tar.xz.part-01 (7 MB, 下载次数 Times of downloads: 11)

2025-7-6 02:13 上传 Uploaded

点击下载
Click to download
Supplementary_Information.tar.xz.part-02 (7 MB, 下载次数 Times of downloads: 11)

2025-7-6 02:14 上传 Uploaded

点击下载
Click to download
Supplementary_Information.tar.xz.part-03 (7 MB, 下载次数 Times of downloads: 11)

2025-7-6 02:15 上传 Uploaded

点击下载
Click to download
Supplementary_Information.tar.xz.part-04 (7 MB, 下载次数 Times of downloads: 10)

2025-7-6 02:16 上传 Uploaded

点击下载
Click to download
Supplementary_Information.tar.xz.part-05 (6.5 MB, 下载次数 Times of downloads: 10)

2025-7-6 02:17 上传 Uploaded

点击下载
Click to download
查看该补充材料的方法稍微有点曲折：
1. 下载完整六个part于同一文件夹，并转移至Linux系统；
2~3. 依次执行下面两条命令——

1. cat Supplementary_Information.tar.xz.part-* > Supplementary_Information.tar.xz
   
2. sha512sum Supplementary_Information.tar.xz

复制代码

——如果返回值为下述字符串则代表成功；

1. 24cf4d8a4ae02a0560c729be03309fdca187cda0b14ad134eef59e986ebb3cf3aff8a97f36c6a308aa4251fbc9f7f71ceb35c825c464c146e9c914b0306cc87b

复制代码

4. 用tar -xvaf Supplementary_Information.tar.xz解压，获得单个.pptx文件；
5. 将.pptx文件转移到Windows系统并用PowerPoint打开，应当看到共7页幻灯片并设置了只读；
6. 点击“仍然编辑”，就可以在后面几页中与3D模型交互了。

当然我理解这么复杂的文件分享流程很容易出问题，如果无法如我所说工作，还请麻烦社长从服务器删除上述附件。

Uus/pMeC6H4-/キ

今天（7月8日）在一楼更新了第2版的脚本，主要是修改了匹配Gaussian输出文本的逻辑以及添加了些提示，请及时更新。总结下目前的适用性测试：
1. 标题所写的IEFPCM和SMD两种最常用的溶剂模型肯定可用，其他较少用且不推荐用的未测试；
2. 最常用的vdW孔洞构建方法肯定可用，二楼幻灯片里展示的较少见的SAS和SES对于Gaussian内置溶剂也可用、对自定义溶剂且溶剂分子半径未知的情形未测试；
3. 有些原子较分散的团簇体系可能构建出不止一个连续的溶质孔洞，此时输出文件会写成

1. Separa: There are  2 disjoint cavities,
   
2. GePol: Number of generator spheres                  =      14
   
3. ...

复制代码

这会干扰第1版的脚本对Number of generator spheres一句的定位，而更新的第2版脚本不再以此句判断溶质孔洞信息结束的位置，因而也可用了；

4. 体系除常规原子外还有Bq鬼原子时，vdW构建的溶质孔洞在该原子处半径为0，例如

1. Spheres list:
   
2. ISph  on   Nord     Re0    Alpha      Xe            Ye            Ze
   
3. ...
   
4.    14  Bq     14    0.0000  1.100      6.000000      0.000000      0.000000

复制代码

第2版脚本会在此提醒一句"Note: sphere 13 has a radius of zero"（脚本内部的列表从0开始计数），最终VMD会显示有个X原子而不存在相应孔洞或顶点，视为可用；

5. 体系除根据原子构建的孔洞外还有自定义添加的孔洞时，比如额外输入写成

1. GeomView
   
2. PrintSpheres
   
3. ExtraSph=1
   
4. {empty}
   
5. 5.0 0.0 0.0 1.0 1.0
   
6. {empty}
   
7. {empty}

复制代码

则产生的自定义孔洞缺少on和Nord两列内容

1. Spheres list:
   
2. ISph  on   Nord     Re0    Alpha      Xe            Ye            Ze
   
3. ...
   
4.    14               1.0000  1.000      5.000000      0.000000      0.000000

复制代码

第2版脚本会提醒"Note: sphere 13 will not be shown due to unrecognized format"，最终VMD也显示不出该孔洞，但GeomView记录的该孔洞表面的顶点仍能显示。

编辑：此外，目前测试的任务均为基态单点能计算，对几何优化或IRC等多步任务solcav.vmd预计会把所有帧的结构与孔洞全合并到一起，而对激发态涉及平衡溶剂/非平衡溶剂和外迭代什么的就难说了，欢迎讨论。