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

Stardust0831 · 发表于 Post on 2025-5-7 03:23:42

本帖最后由 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盘。）

脚本在此下载：

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

ESP_orca.sh (5.62 KB, 下载次数 Times of downloads: 139)

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

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*物理内存/使用核数。

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

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

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所在目录加入了系统环境变量，可以直接输入：

vmd.exe

复制代码

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

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刻度的上下限，比如在命令行中输入这句：

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的图片。

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

本文转载请注明出处。

yubozou@163.com · 发表于 Post on yesterday 15:55

请问下我在计算时出现以下报错是什么原因，我需要怎么调整？λ bash ESP_orca.sh
Geting the structure from PANIPA.in ...
~~~~~~~~~~~~We are runing pre-optimization~~~~~~~~~~~~
Running PANIPA_preopt.inp ...
TOTAL RUN TIME: 0 days 0 hours 2 minutes 44 seconds 734 msec

Geting the structure from PANIPA_preopt.out ...
~~~~~~~~~~~~We are runing optimization~~~~~~~~~~~~
Running PANIPA_opt.inp ...
TOTAL RUN TIME: 1 days 4 hours 10 minutes 54 seconds 993 msec

Geting the structure from PANIPA_opt.out ...
~~~~~~~~~~~~We are runing single point~~~~~~~~~~~~
Running PANIPA_sp.inp ...

[file orca_tools/SHARK/shark_setup_ri.cpp, line 487, Process 2]:

WARNING! Potentially linear dependencies in the auxiliary basis
===>> At the moment, SHARK is not equipped to deal with this
... sorry, have to bail out
[file orca_tools/SHARK/shark_setup_ri.cpp, line 487, Process 0]:

WARNING! Potentially linear dependencies in the auxiliary basis
===>> At the moment, SHARK is not equipped to deal with this
... sorry, have to bail out
[file orca_tools/SHARK/shark_setup_ri.cpp, line 487, Process 3]:

WARNING! Potentially linear dependencies in the auxiliary basis
===>> At the moment, SHARK is not equipped to deal with this
... sorry, have to bail out
[file orca_tools/SHARK/shark_setup_ri.cpp, line 487, Process 1]:

WARNING! Potentially linear dependencies in the auxiliary basis
===>> At the moment, SHARK is not equipped to deal with this
... sorry, have to bail out
[file orca_tools/qcmsg.cpp, line 394]:
  .... aborting the run

Reading the input file PANIPA_sp.gbw             ... done
Writing the output file PANIPA_sp.molden.input    ... done
~~~~~~~~~~~~We are analyzing the molden file~~~~~~~~~~~~
Export cub file from PANIPA_sp.molden.input ...
forrtl: severe (24): end-of-file during read, unit 10, file D:\software\esp\test\ANFPANI\PANIPA_sp.molden.input
Image             PC             Routine          Line       Source
Multiwfn.exe    00007FF72742B7DF  Unknown             Unknown  Unknown
Multiwfn.exe    00007FF72740C8B3  Unknown             Unknown  Unknown
Multiwfn.exe    00007FF726FDCDED  Unknown             Unknown  Unknown
Multiwfn.exe    00007FF726F90F20  Unknown             Unknown  Unknown
Multiwfn.exe    00007FF7272B0DE0  Unknown             Unknown  Unknown
Multiwfn.exe    00007FF727DA274E  Unknown             Unknown  Unknown
Multiwfn.exe    00007FF72806F8D4  Unknown             Unknown  Unknown
KERNEL32.DLL    00007FF81645259D  Unknown             Unknown  Unknown
ntdll.dll       00007FF816DEAF78  Unknown             Unknown  Unknown
mv: cannot stat 'density.cub': No such file or directory
mv: cannot stat 'totesp.cub': No such file or directory

Analysing PANIPA_sp.molden.input ...
forrtl: severe (24): end-of-file during read, unit 10, file D:\software\esp\test\ANFPANI\PANIPA_sp.molden.input
Image             PC             Routine          Line       Source
Multiwfn.exe    00007FF72742B7DF  Unknown             Unknown  Unknown
Multiwfn.exe    00007FF72740C8B3  Unknown             Unknown  Unknown
Multiwfn.exe    00007FF726FDCDED  Unknown             Unknown  Unknown
Multiwfn.exe    00007FF726F90F20  Unknown             Unknown  Unknown
Multiwfn.exe    00007FF7272B0DE0  Unknown             Unknown  Unknown
Multiwfn.exe    00007FF727DA274E  Unknown             Unknown  Unknown
Multiwfn.exe    00007FF72806F8D4  Unknown             Unknown  Unknown
KERNEL32.DLL    00007FF81645259D  Unknown             Unknown  Unknown
ntdll.dll       00007FF816DEAF78  Unknown             Unknown  Unknown
mv: cannot stat 'surfanalysis.pdb': No such file or directory

mv: cannot stat '*.cub': No such file or directory
mv: cannot stat '*_surfanalysis.pdb': No such file or directory
Info) VMD for WIN64, version 1.9.3 (May 12, 2021)
Info) http://www.ks.uiuc.edu/Research/vmd/
Info) Email questions and bug reports to vmd@ks.uiuc.edu
Info) Please include this reference in published work using VMD:
Info) Humphrey, W., Dalke, A. and Schulten, K., `VMD - Visual
Info) Molecular Dynamics', J. Molec. Graphics 1996, 14.1, 33-38.
Info) -------------------------------------------------------------
Info) Multithreading available, 8 CPUs detected.
Info) Free system memory: 10GB (63%)
Info) No CUDA accelerator devices available.
Info) OpenGL renderer: Intel(R) Iris(R) Xe Graphics
Info) Features: STENCIL MDE CVA MTX NPOT PP PS GLSL(OVFGS)
Info) Full GLSL rendering mode is available.
Info) Textures: 2-D (16384x16384), 3-D (1024x1024x1024), Multitexture (8)
Info) No joysticks found.  Joystick interface disabled.
couldn't read file "showorb.vmd": no such file or directory
-10
10
no files matched glob pattern "*density.cub"
1
can't read "allfiles": no such variable
molinfo: get: no molecule exists with id -1

yubozou@163.com · 发表于 Post on 2026-2-23 10:35:34

Stardust0831 发表于 2026-2-22 02:11
tachyon_WIN64.exe是在vmd安装文件夹下的，需要确保1.使用的确实是64位1.9.3的vmd，而非官网上那个32位的 ...

感谢感谢，原来我的是32位的

Stardust0831 · 发表于 Post on 2026-2-22 02:11:30

yubozou@163.com 发表于 2026-2-18 16:53
你好，我用上面新的脚本运行后vmd窗口正常，但是.bmp格式的文件没办法输出，窗口提示couldn't execute "tac ...

tachyon_WIN64.exe是在vmd安装文件夹下的，需要确保1.使用的确实是64位1.9.3的vmd，而非官网上那个32位的版本2.已经把vmd所在目录加到了系统环境变量里

yubozou@163.com · 发表于 Post on 2026-2-18 16:53:12

你好，我用上面新的脚本运行后vmd窗口正常，但是.bmp格式的文件没办法输出，窗口提示couldn't execute "tachyon_WIN64.exe": no such file or directory，想问下要怎么操作，是我哪里有问题？谢谢

CM234 · 发表于 Post on 2026-2-12 18:45:49

Stardust0831 发表于 2026-1-29 16:51
当前文件夹下有哪些文件呢？

感谢回复，已经排除了，我有一个路径写错了。

Stardust0831 · 发表于 Post on 2026-1-29 16:51:35

CM234 发表于 2026-1-23 17:32
请问我将vmd，orca和Multiwfn都加入系统变量，也修改了ESP_orca.sh文件中的路径，运行后出现如下报错，请问 ...

当前文件夹下有哪些文件呢？

CM234 · 发表于 Post on 2026-1-23 17:32:45

本帖最后由 CM234 于 2026-1-23 18:22 编辑

请问我将vmd，orca和Multiwfn都加入系统变量，也修改了ESP_orca.sh文件中的路径，运行后出现如下报错，请问要怎么排除呢？附上我的输入文件

wwwwwt · 发表于 Post on 2025-8-15 08:39:00

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

感谢感谢！

wwwwwt · 发表于 Post on 2025-8-15 08:38:51

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

感谢回复！

Stardust0831 · 发表于 Post on 2025-8-14 21:29:28

本帖最后由 Stardust0831 于 2025-9-24 13:53 编辑

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//.fchk/_density.cub}
mv -f totesp.cub ${inf//.fchk/_ESP.cub}
echo
done

复制代码

4.绘制esp，此时基本上类似《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图（含视频）》里的流程，只是.vmd文件要换成这个：

ESP_stardust0831.vmd (2.53 KB, 下载次数 Times of downloads: 6)
和原版的差异是cub文件命名是按“分子名_ESP.cub”这种格式的，可以直接看cub认出是哪个分子，并且有我在二楼提到的一些自定义函数，比如“xr 文件名”直接按特定文件名导出并渲染图片等。

Stardust0831 · 发表于 Post on 2025-8-14 21:07:49

本帖最后由 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就可以读取了。

wwwwwt · 发表于 Post on 2025-8-14 09:47:05

楼主你好，我在运行esp_sh时，一直出现输入文件无效的错误，附件是我的输入文件，请问是哪里出错了么？

wwwwwt · 发表于 Post on 2025-8-14 09:43:38

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

感谢

wal · 发表于 Post on 2025-8-13 14:59:54

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

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

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[辅助/分析程序] 使用Multiwfn+ORCA自动批量计算静电势与偶极矩的脚本

评分 Rate

评分 Rate