gau_chemsh.py: Tcl-ChemShell调用Gaussian优化器
gau_chemsh.py: Enabling Tcl-ChemShell Using Gaussian Optimizer
Jiawei Xu
Released: 2024-05-09 / Updated: 2024-05-09
Tcl-ChemShell（以下简称chemsh）是常用的QM/MM计算接口程序。默认情况下，chemsh使用DL-find作为优化器。chemsh支持做极小点优化、使用dimer或NEB (nudged elastic band)方法做过渡态搜索，以及最低能量路径(minimum energy path, MEP)搜索。实践证明，虽然一般认为dimer方法做过渡态搜索是很高效的，但其结果对初始结构要求比较高，对于复杂体系一般需要通过柔性扫描给出起始的两个dimer结构，而失败概率仍比较高。另外，chemsh不支持IRC计算。本文将介绍gau_chemsh.py的使用方法。
1. 柔性扫描
chemsh有自带的surface关键词做柔性扫描，但没有简便的方法可以生成柔性扫描所需的冗余内坐标矩阵(z-matrix)，也不方便添加被扫描的变量。虽然可以在输入文件中编写复杂的tcl脚本用于更新结构，但这只适用于处理被扫描变量是原子距离的情况，无法处理更复杂的结构变量。本例我们介绍使用gau_chemsh.py为chemsh调用Gaussian优化器做柔性扫描计算的方法。所用软件版本为Tcl-ChemShell 3.7.1和Gaussian16 Revision C.01。
1.1 输入文件
我们使用chemsh作为ORCA和DL_poly的接口，其中参数使用AMBER格式。准备chemsh做QM/MM计算所必须的输入文件scan.c、scan.prmtop和scan.chm，其中scan.c为起始结构，scan.prmtop为整个体系的参数文件；scan.chm为chemsh输入文件的模板，但不包含具体计算命令，其具体内容如下：
set sysname scan
set program "~/programs/orca-5.0.4/orca"
source ~/files/tools/orca5.0-chemsh.tcl
set qm_atom_numbers { 1-33 }
set active_atom_numbers { 1-61 }
set qm_charge -2
set qm_mult 1
set orcasimpleinput "! xtb"
set orcablocks { %pal nprocs 16 end }
set qmmm [ list \
coupling = shift \
qm_region = ${qm_atom_numbers} \
qm_theory = orca : [ list \
executable = ${program} \
jobname = ${sysname}_orca \
orcasimpleinput = ${orcasimpleinput} \
orcablocks = ${orcablocks} \
charge = ${qm_charge} \
mult = ${qm_mult} ] \
mm_theory = dl_poly : [ list \
amber_prmtop_file = ${sysname}.prmtop \
exact_srf = yes \
use_pairlist = no \
mxlist = 40000 \
mxexcl = 2000 \
cutoff = 1000 \
debug_memory = no \
scale14 = [ list [ expr 1 / 1.2 ] 0.5 ] \
save_dl_poly_files = yes \
conn = ${sysname}.c \
list_option = none ]]上述内容包含了体系的相关设置，以及变量名为qmmm的能量梯度方法设定，这将被gau_chemsh.py用于在各种计算任务中设定“theory = hybrid : ${qmmm}”。注意，本模板中使用“qm_atom_numbers”和“active_atom_numbers”设定QM区域和活性区域的原子序号，这两个变量名不能改变，gau_chemsh.py将据此确定QM区域和活性区域的范围。
另外准备gau_chemsh.py的输入文件scan.inp，其具体内容如下：
# input file for gau_chemsh.py
c0 scan.c # necessary: input structure
chm scan.chm # necessary: chemsh input file template
sysname scan # necessary: point to ${sysname}.e, ${sysname}.g and ${sysname}.h上述输入文件包含了三个关键词：(1) c0：指定初始结构；(2) chm：指定chemsh输入模板文件；(3) sysname：设置体系名称，这将用于一系列文件的命名。为了使用户使用更为灵活，gau_chemsh.py不会从chemsh模板文件中尝试读取初始结构或体系名称，因此以上关键词都是必须设定的。
1.2 柔性扫描计算
在命令行中运行：python gau_chemsh.py -i scan.inp –init，生成用于Gaussian做柔性扫描任务的输入文件scan.gjf。默认产生的是做过渡态搜索计算的输入文件，手动修改相应关键词，即指定opt=modredundant，并添加被扫描变量的冗余内坐标设定。对于本例，其具体内容如下：
! Tcl-ChemShell TS Task Using Gaussian16 Optimizer
%nprocshared=1
#p opt(modredundant,nomicro) nosymm ugbs
external="python gau_chemsh.py -i scan.inp --main"
generated by XMECP
0 1
O -2.2513047912 -2.6780736895 -0.4004817984
O -1.3880971945 -0.5198143980 -0.5873724977
... ...
O 1.2902934949 3.6770722855 7.7870162694
B 2 9 S 10 0.1
可以根据实际需要增减关键词。注意：如果添加了calcfc、recalc=n或calcall，gau_chemsh.py会按需进行二阶导数计算，但chemsh不支持静电嵌入的QM/MM解析二阶导数，因此是计算数值二阶导数，耗时很高，gau_chemsh.py默认只对QM区域原子计算数值二阶导数。默认的关键词中，nomicro和nosymm不能删去，ugbs不建议删去。另外，不要修改并行核数，否则Gaussian会总是占用一半的核数而不执行任何计算。
如果当前系统中使用“python”命令所调用的python不适用于gau_chemsh.py，例如版本过低（需要python3）或无相关库支持等，可以自行修改external中运行python的命令，使之对应于合适的python可执行文件。
完成以上工作后，运行：g16 scan.gjf，即进行柔性扫描计算。对于集群计算的情况，则是提交一个Gaussian作业，并注意提交队列系统的脚本中应包含所有用到软件的环境变量，即在单独做相同chemsh作业的基础上还额外需要Gaussian的环境变量。计算完成后，可使用GaussView直接查看柔性扫描的结果。该方法适用于扫描键长、键角和二面角，也适用于做各种限制性优化。对于Gaussian16 A.03及以后版本，还可以通过定义广义内坐标(generalized internal coordinate, GIC)来引入更为复杂的结构变量，例如两个芳香环之间的距离、两个分子质心之间的距离等。