计算化学公社

标题: 使用Multiwfn+ORCA自动批量计算静电势与偶极矩的脚本 [打印本页]

作者
Author: Stardust0831 时间: 2025-5-7 03:23
标题: 使用Multiwfn+ORCA自动批量计算静电势与偶极矩的脚本
本帖最后由 Stardust0831 于 2025-5-9 20:08 编辑

静电势（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盘。）

脚本在此下载： (, 下载次数 Times of downloads: 6) (已过时)
5月9日更新： (, 下载次数 Times of downloads: 49)

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

export Multiwfnpath=/g/Downloads/Multiwfn_3.8_dev_bin_Win64-20250220/Multiwfn_3.8_dev_bin_Win64
export PATH=$PATH:/g/Downloads/Multiwfn_3.8_dev_bin_Win64-20250220/Multiwfn_3.8_dev_bin_Win64
orca='/d/orca6.0.1/orca.exe'
File_extension="gjf"
nprocs=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*物理内存/使用核数。

本脚本的执行效果如下：
新建一个文件夹，将结构文件和脚本放入其中：
(, 下载次数 Times of downloads: 15)
打开cmder，cd到对应目录，然后输入：

bash ESP_orca.sh

复制代码

(, 下载次数 Times of downloads: 12)
可以看到脚本自动调用的multiwfn生成了我们需要的输入文件并调用orca执行了计算半经验级别的预优化、dft级别的优化、高精度的单点，并调用multiwfn进行了后处理。计算结束后，自动把最重要的完成了文件的归类：ESP文件夹是绘制ESP需要用的文件，opt文件夹是结构优化的输出文件（建议计算完成后再检查一下是否有虚频），sp是单点计算的输出文件，tmp是一些我认为不常用的、可以删了的文件（以防你可能需要用到这些文件，我选择放入特定文件夹而非直接删）。
当前目录下还产生了surface_analysis.txt和DIPOLE_MOMENT.txt（更新后的脚本会把偶极矩和ESP上下限合并输出到Information.txt这个文件里）文件文件，分别记录了每个体系的表面静电势的统计情况和偶极矩：
(, 下载次数 Times of downloads: 15)
每段数据前会先输出文件名。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 文件名”的方式即可完成导出与渲染，见一楼回帖）
(, 下载次数 Times of downloads: 15)
在ESP文件夹下右键打开powershell（或者打开powershell以后自行cd过去）：
(, 下载次数 Times of downloads: 15) (, 下载次数 Times of downloads: 15)
输入vmd.exe所在路径，如果将vmd所在目录加入了系统环境变量，可以直接输入：

vmd.exe

复制代码

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

source ESP_stardust0831.vmd

复制代码

可以加载当前目录下的所有cub文件
(, 下载次数 Times of downloads: 12)
在vmd的main窗口中，有T A D F四个字母，T表示Top，双击对应的位置可以选择处理的模型。D代表显示、F代表冻结（均可以双击来改变状态）
(, 下载次数 Times of downloads: 16)
我们先处理H2O的模型，把T点到H2O模型的前面，然后再命令行输入top三个字母并回车，可以实现只显示被Top的模型并调整视角到合适的位置。（top是我自定义的命令）
(, 下载次数 Times of downloads: 15)
然后加上色彩刻度轴，参考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刻度的上下限，比如在命令行中输入这句：

mol scaleminmax top 1 -30 35

复制代码

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

render Tachyon vmdscene.dat

复制代码

把当前的可视化效果导出成vmdscene.dat文件，双击Render_stardust0831.bat就可以把数据从vmdscene.dat中导出成vmdscene.bmp的图片。
(, 下载次数 Times of downloads: 13)
把vmdscene.bmp文件手动改成体系的名字（防止下一次画图的时候被覆盖）

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

draw color green
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，算是给我的最大鼓励。有任何问题可以留言，我会常看这个帖子的。

本文转载请注明出处。

作者
Author: Stardust0831 时间: 2025-5-7 19:47
本帖最后由 Stardust0831 于 2025-5-7 21:28 编辑

我自己占一下沙发，方便后续放更新。
脚本正在修改，目前已实现在结构优化后发现有虚频时将虚频输出
(, 下载次数 Times of downloads: 14)

更新已完成。 (, 下载次数 Times of downloads: 6)
美化了输出结果同时如果出现虚频会将虚频输出到屏幕上。将ESP最大和最小值、偶极矩大小和方向输出到了Information.txt，使数据更紧凑。
(, 下载次数 Times of downloads: 15) (, 下载次数 Times of downloads: 14)
计算完成后会自动调用vmd并加载绘图脚本。
增加了以下自定义函数以方便绘图：
top命令：让vmd只显示top体系并调整到合适初始视角。（载入所有文件后后默认会加载一次top命令）

proc top {} {
mol off all
mol on top
display resetview
}

复制代码

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

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

复制代码

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

xr H2O

复制代码

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

proc esp {colorlow colorhigh} {
mol scaleminmax top 1 $colorlow $colorhigh
}

复制代码

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

esp -15 24

复制代码

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

ou 0.87365130 0.50440282 -0.35666665 3

复制代码

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

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

作者
Author: wwwwwt 时间: 前天 11:08
请问有针对Gaussian16 Linux版本的批量计算静电势脚本吗？

作者
Author: wwwwwt 时间: 前天 14:31
请问显示静电势极大值和极小值时，是否能显示其极大极小值在哪个原子或基团上么？

作者
Author: wal 时间: 前天 14:59

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

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

作者
Author: wwwwwt 时间: yesterday 09:43

wal 发表于 2025-8-13 14:59
显示在哪个原子/基团上有点难吧。顶多告诉你index是多少，具体属于哪个原子/基团得自己判断

感谢

作者
Author: wwwwwt 时间: yesterday 09:47
楼主你好，我在运行esp_sh时，一直出现输入文件无效的错误，附件是我的输入文件，请问是哪里出错了么？

作者
Author: Stardust0831 时间: yesterday 21:07
本帖最后由 Stardust0831 于 2025-8-14 21:11 编辑

wwwwwt 发表于 2025-8-14 09:47
楼主你好，我在运行esp_sh时，一直出现输入文件无效的错误，附件是我的输入文件，请问是哪里出错了么？

gjf是gaussian的格式，我会用gjf当输入是因为我建模的主力软件是gview。
你文件里写的事实上是orca的格式和关键词（并且关键词是过时的），这显然和gjf这个后缀不符，multiwfn当gaussian输入文件来读取自然认不出来。
把脚本里的

File_extension="gjf"

复制代码

改成

File_extension="in"

复制代码

并且把这个orca格式的文件的后缀改成.in，这样multiwfn就可以读取了。

作者
Author: Stardust0831 时间: yesterday 21:29
本帖最后由 Stardust0831 于 2025-8-14 21:32 编辑

wwwwwt 发表于 2025-8-13 11:08
请问有针对Gaussian16 Linux版本的批量计算静电势脚本吗？

g16用不了GFN2-xTB这种半经验级别，所以我自己的工作流不会用g16。
简单给你讲一下怎么批量运行：
（假设gaussian已经被正确安装，在环境变量里了）
1.批量执行gjf文件

for inf in *.gjf
do
echo Running ${inf} ...
g16 < ${inf} > ${inf//gjf/out}
done

复制代码

这5行单独当一个sh脚本，可以执行当前目录下所有以.gjf结尾的文件。

2.批量生成波函数文件：

for inf in *.chk
do
formchk ${inf}
done

复制代码

3.生成用于绘图的格点文件：

for inf in *.fchk
do
echo Export cub file from ${inf} ...
Multiwfn.exe $inf << EOF > log.txt
5
1
3
2
0
5
12
1
2
0
r
totesp.cub
13
11
5
627.5095
0
totesp.cub
-1
q
EOF
mv -f density.cub ${inf//.molden.input/_density.cub}
mv -f totesp.cub ${inf//.molden.input/_ESP.cub}
echo
done

复制代码

4.绘制esp，此时基本上类似《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图（含视频）》里的流程，只是.vmd文件要换成这个：
(, 下载次数 Times of downloads: 1)
和原版的差异是cub文件命名是按“分子名_ESP.cub”这种格式的，可以直接看cub认出是哪个分子，并且有我在二楼提到的一些自定义函数，比如“xr 文件名”直接按特定文件名导出并渲染图片等。

作者
Author: wwwwwt 时间: 9 hour ago

Stardust0831 发表于 2025-8-14 21:07
gjf是gaussian的格式，我会用gjf当输入是因为我建模的主力软件是gview。
你文件里写的事实上是orca的格 ...

感谢回复！

作者
Author: wwwwwt 时间: 9 hour ago

Stardust0831 发表于 2025-8-14 21:29
g16用不了GFN2-xTB这种半经验级别，所以我自己的工作流不会用g16。
简单给你讲一下怎么批量运行：
（假 ...

感谢感谢！

欢迎光临计算化学公社 (http://bbs.keinsci.com/)