使用Multiwfn+ORCA自动批量计算静电势与偶极矩的脚本

Stardust0831

静电势（ESP）和偶极矩作为非常常用的数据，批量计算的需求很旺盛。
其中，静电势对于考察分子间静电相互作用、预测反应位点、预测分子性质等方面有重要意义，相关学习资料见《静电势与平均局部离子化能综述合集》

本脚本旨在自动调用Multiwfn和ORCA自动实现利用刚建模完成的结构文件（如高斯的gjf文件，其实任何multiwfn可以识别的结构文件都可以），自动调用ORCA完成GFN2-xTB级别的预优化，B3LYP-D3(BJ)/def2-SVP级别的结构优化和B3LYP-D3(BJ)/def2-TZVPD级别的单点；并调用multiwfn自动从最后单点得到的波函数文件计算得到静电势的格点文件、以及从单点的输出文件中提取偶极矩信息。
有别于之前的脚本，此处生成的格点文件的文件名会包含体系名，方便辨认是哪个体系的；计算完成后也会自动归档文件；并设定了vmd的快速导出并渲染的自定义指令。

本脚本编写过程中参考了以下博文/帖子：《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图（含视频）》《使用Multiwfn结合VMD分析和绘制分子表面静电势分布》《使用Multiwfn计算分子片段的偶极矩和复合物中单体的偶极矩》
本脚本涉及tcl编写的部分是通过北京科音分子动力学与GROMACS培训班学习的。

本脚本需要使用的软件有ORCA6，Multiwfn和VMD1.9.3（64位），安装教程可参考：
ORCA：《量子化学程序ORCA的安装方法》（请使用最新版，如是旧版本、需解决otool_xtb调用问题、或自行去除脚本中关于GFN2-xTB级别的预优化相关部分）（windows用户需参考此博文完成cmder安装并设置bash）“利用cmder可以令ORCA在Windows下的使用明显更方便，cmder是一个第三方的文本终端。首先去https://cmder.app下载Full版的cmder，然后解压到你平时安装应用程序的目录。启动cmder窗口后，在cmder的标题位置点右键选Settings，在General页面里选择{bash::bash as Admin}，然后把cmder关了。从此之后，新开的cmder终端里的命令写法就和Linux的Bash环境下一样了，连awk、vi等常用工具都有。”
Multiwfn：http://sobereva.com/multiwfn/download.html（windows端下载解压即用，linux系统参考Multiwfn在Linux下安装的中文说明）
VMD：使用此处提供的安装包安装vmd1.9.3的64位版，链接: https://pan.baidu.com/s/1J0hfvqjrVYdDM_iW_eqagQ?pwd=e229 提取码: e229（推荐安装完成后将vmd目录加入系统环境变量，请尽量避免将vmd安装在C盘。）

脚本在此下载： ESP_orca.sh (5.05 KB, 下载次数 Times of downloads: 6)

2025-5-7 03:11 上传 Uploaded

点击下载
Click to download
阅读权限: 10

(已过时)
5月9日更新： ESP_orca.sh (5.62 KB, 下载次数 Times of downloads: 44)

2025-5-9 19:52 上传 Uploaded

点击下载
Click to download

使用前需要修改脚本开头的环境变量：

1. export Multiwfnpath=/g/Downloads/Multiwfn_3.8_dev_bin_Win64-20250220/Multiwfn_3.8_dev_bin_Win64
   
2. export PATH=$PATH:/g/Downloads/Multiwfn_3.8_dev_bin_Win64-20250220/Multiwfn_3.8_dev_bin_Win64
   
3. orca='/d/orca6.0.1/orca.exe'
   
4. File_extension="gjf"
   
5. nprocs=6
   
6. maxcore=8000

复制代码

把/g/Downloads/Multiwfn_3.8_dev_bin_Win64-20250220/Multiwfn_3.8_dev_bin_Win64/Multiwfn.exe换成你的multiwfn路径
把/d/orca6.0.1/orca.exe换成你的orca.exe所在位置
File_extension="gjf"表示使用的结构文件是以gjf为后缀的。
nprocs=6后面的数值6表示使用的核数，改成你的。这个值不可以超过cpu的物理核心数。
maxcore=8000表示每个核心用8000MB内存，不可以超过0.8*物理内存/使用核数。

本脚本的执行效果如下：
新建一个文件夹，将结构文件和脚本放入其中：

打开cmder，cd到对应目录，然后输入：

1. bash ESP_orca.sh

复制代码

可以看到脚本自动调用的multiwfn生成了我们需要的输入文件并调用orca执行了计算半经验级别的预优化、dft级别的优化、高精度的单点，并调用multiwfn进行了后处理。计算结束后，自动把最重要的完成了文件的归类：ESP文件夹是绘制ESP需要用的文件，opt文件夹是结构优化的输出文件（建议计算完成后再检查一下是否有虚频），sp是单点计算的输出文件，tmp是一些我认为不常用的、可以删了的文件（以防你可能需要用到这些文件，我选择放入特定文件夹而非直接删）。
当前目录下还产生了surface_analysis.txt和DIPOLE_MOMENT.txt（更新后的脚本会把偶极矩和ESP上下限合并输出到Information.txt这个文件里）文件文件，分别记录了每个体系的表面静电势的统计情况和偶极矩：

每段数据前会先输出文件名。esp的最大值和最小值可直接在surface_analysis.txt查看Minimal value: 和Maximal value: 得到；DIPOLE_MOMENT.txt中的Magnitude (Debye)  就是偶极矩的值，单位是debye， Total Dipole Moment    :   是偶极矩的三个分量（绘制偶极矩箭头需要用的）。

ESP文件夹中产生了_density.cub和_ESP.cub分别是电子密度的格点和静电势的格点，画ESP图要用，_surfanalysis.pdb是ESP的极值点。此外给的ESP_stardust0831.vmd是vmd的绘图脚本，Render_stardust0831.bat是图片渲染脚本（给vmd目录设了环境变量以后可以直接双击渲染在vmd中绘制导出的dat文件）。这些文件会由脚本自动生成。
（更新后的脚本不再生成bat渲染脚本，直接在vmd中用“xr 文件名”的方式即可完成导出与渲染，见一楼回帖）

在ESP文件夹下右键打开powershell（或者打开powershell以后自行cd过去）：

输入vmd.exe所在路径，如果将vmd所在目录加入了系统环境变量，可以直接输入：

1. vmd.exe

复制代码

然后就可以进入vmd了。在vmd中输入：

1. source ESP_stardust0831.vmd

复制代码

可以加载当前目录下的所有cub文件

在vmd的main窗口中，有T A D F四个字母，T表示Top，双击对应的位置可以选择处理的模型。D代表显示、F代表冻结（均可以双击来改变状态）

我们先处理H2O的模型，把T点到H2O模型的前面，然后再命令行输入top三个字母并回车，可以实现只显示被Top的模型并调整视角到合适的位置。（top是我自定义的命令）

然后加上色彩刻度轴，参考196博文中的这个方法

之后给图上加上色彩刻度轴。选Extensions-Visualization-Color Scale Bar，Color bar width设为0.08，Display title选on并且将Color bar title里写上ESP (kcal/mol)，Minimum和Maximum scale value分别填-22和22，Number of axis labels输入10，Color labels选Black，Label format选Decimal。然后点Draw Color Scale Bar按钮，色彩刻度就出现在画面中了，并且VMD Main窗口中多出了一个名为Color Scale Bar的一项。然后调整它的大小和位置，即双击VMD Main窗口中Color Scale Bar那一项当中的F标签使之变为红色（即不让色彩刻度轴在画面中的位置冻结），而双击其它项目的F标签使它们的F变为黑色（让它们的位置冻结住）。然后点击VMD的OpenGL图形窗口激活之，按t键进入平移模式，然后拖动鼠标将色彩刻度轴放置到合适位置，并且用鼠标滚轮调整它的大小。调合适之后再按r键恢复旋转视角模式，并且在VMD Main里将Color Scale Bar那一项的F重新双击成黑色，而其它三项的F重新双击为红色。

应当注意我此处把Autoscale也开了，使上下限的刻度符合我的调整。

简单调整模型的ESP刻度的上下限，比如在命令行中输入这句：

1. mol scaleminmax top 1 -30 35

复制代码

这个就表示把当前的ESP的刻度下限调整到-30kcal/mol，上限调整到35kcal/mol。
调整妥当后，在vmd中输入

1. render Tachyon vmdscene.dat

复制代码

把当前的可视化效果导出成vmdscene.dat文件，双击Render_stardust0831.bat就可以把数据从vmdscene.dat中导出成vmdscene.bmp的图片。

把vmdscene.bmp文件手动改成体系的名字（防止下一次画图的时候被覆盖）

也可以使用使用Multiwfn计算分子片段的偶极矩和复合物中单体的偶极矩中的类似这种之类指令绘制偶极矩，1.227306 -0.128087 0.650833 是偶极矩方向，最后的2是系数、系数越大箭头越长。（后来我对这一功能进行了一定修改，见一楼回帖）

1. draw color green
   
2. drawarrow all 1.227306 -0.128087 0.650833 2

复制代码

文末附上文中演示时用到的文件：
链接: https://pan.baidu.com/s/1p0hpTVwDzdu9PWz8_5raRw?pwd=96wp 提取码: 96wp

自动调用orca和multiwfn部分的脚本是我纯手写的。绘制ESP部分的代码是在社长的基础上做了修改，也加入了额外的自定义命令来方便批量绘制的这一过程。为了偶极矩箭头的绘图方便，我将偶极矩箭头功能是作为附加功能直接加了进去，偶极矩箭头的绘制的这部分部分的代码Written by Tian Lu。

脚本好用的话，可以给我加点eV，算是给我的最大鼓励。有任何问题可以留言，我会常看这个帖子的。

本文转载请注明出处。

Stardust0831

我自己占一下沙发，方便后续放更新。
脚本正在修改，目前已实现在结构优化后发现有虚频时将虚频输出



更新已完成。 ESP_orca.sh (5.62 KB, 下载次数 Times of downloads: 4)

2025-5-7 21:27 上传 Uploaded

点击下载
Click to download

美化了输出结果同时如果出现虚频会将虚频输出到屏幕上。将ESP最大和最小值、偶极矩大小和方向输出到了Information.txt，使数据更紧凑。

计算完成后会自动调用vmd并加载绘图脚本。
增加了以下自定义函数以方便绘图：
top命令：让vmd只显示top体系并调整到合适初始视角。（载入所有文件后后默认会加载一次top命令）

1. proc top {} {
   
2. mol off all
   
3. mol on top
   
4. display resetview
   
5. }

复制代码

xr命令：直接导出并渲染。

1. proc xr {file_name} {
   
2. render Tachyon $file_name.dat
   
3. exec tachyon_WIN64.exe $file_name.dat -aasamples 24 -trans_vmd -res 1024 768 -format BMP -o $file_name.bmp
   
4. }

复制代码

比如这句命令就是将当前vmd可视化窗口的显示直接导出成H2O.bmp文件

1. xr H2O

复制代码

esp命令，快速调整绘图的上下限：

1. proc esp {colorlow colorhigh} {
   
2. mol scaleminmax top 1 $colorlow $colorhigh
   
3. }

复制代码

比如这句命令就是把ESP的色彩下限改成-15kcal/mol，上限改成24kcal/mol

1. esp -15 24

复制代码

ou命令，参考自《使用Multiwfn计算分子片段的偶极矩和复合物中单体的偶极矩》
使用这种写法可以画出朝 0.87365130    0.50440282   -0.35666665方向（方向从上文提到的Information.txt中获取）箭头长度乘3倍的偶极矩箭头

1. ou  0.87365130    0.50440282   -0.35666665 3

复制代码

turbo命令，参考自《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图（含视频演示）》，如果手动修改脚本使用vmd1.9.4，可以绘制出彩虹色的色彩刻度轴。

KuangL

成功复现！
总的来说配置环境的时间并不长，跟着教程一步一步来一会儿就完成了
确实是一个很方便省事的小脚本，大大简化了绘制ESP的难度
期待下一次更新！

wwwwwt

请问有针对Gaussian16 Linux版本的批量计算静电势脚本吗？

wwwwwt

请问显示静电势极大值和极小值时，是否能显示其极大极小值在哪个原子或基团上么？

wal

wwwwwt 发表于 2025-8-13 14:31
请问显示静电势极大值和极小值时，是否能显示其极大极小值在哪个原子或基团上么？

显示在哪个原子/基团上有点难吧。顶多告诉你index是多少，具体属于哪个原子/基团得自己判断