GKS-EDA计算简介

2022-05-17 14:59:11 浏览数 (1)

GKS-EDA是一种常见的能量分解方法，用于将结合能分解为各个组分，便于比较类似体系间各能量组分分别有多大影响。由于是在DFT水平下做计算，可计算的体系颇大（相比于PSI4里的SAPT2 /aug-cc-pVDZ和SAPT2 (3)δMP2/aug-cc-pVTZ而言）。有不少同学询问笔者如何使用及经常碰到SCF不收敛问题，因此写一篇介绍一下。笔者不是GKS-EDA开发人员，不能保证该文安装过程与将来官方发布的一致。若以后官方更新程序或更新安装方式，不保证此文步骤仍适用。有不妥之处还请各位EDA老司机指正。

1. 安装

GKS-EDA计算由XEDA程序完成，但其实XEDA目前并非独立程序，而是一个修改版的GAMESS程序，所以笔者看官方近年来使用“GAMESS中的XEDA模块”之类的称呼，确实更为合理。以前该程序的获取是向厦门大学苏培峰老师发邮件申请，然后他会发给你一个压缩包，里面是编译好的程序，解压后修改一下环境变量即可使用，没有编译过程。如今获取XEDA程序无需通过邮件往来，可以自己在官网申请获取。先到官网

代码语言：javascript复制
http://xmvb.xmu.edu.cn

注册一个账号，然后到下载页

代码语言：javascript复制
http://xmvb.xmu.edu.cn/downloads

申请获取。笔者申请完了看几分钟内没动静，就过了几天再登录，发现个人中心里的Program Download页面就可以直接下载了。下载下来是一个xeda-patch.tar.gz压缩包，里面是修改代码的脚本xeda-patch、代码文件夹和说明文件README.md。这个意思是，我们需要先去GAMESS官网下载GAMESS源代码，然后使用xeda-patch脚本修改GAMESS源代码，最后编译这个修改版的GAMESS。显然，这与以前通过邮件获取二进制压缩包的步骤不一样。GAMESS官网为

代码语言：javascript复制
https://www.msg.chem.iastate.edu/gamess/License_Agreement.html

有些同学会使用旧版GAMESS，不过笔者比较推荐2021-R2版本的GAMESS，因为最新的GAMESS支持赝势L=5的h角动量，这在算镧系、锕系等最后两周期的金属时经常会用到，若使用低于2021的版本在这个功能上会报错。如果不清楚版本号，可以打开gamess/RELEASE.md文件查看。由于要修改GAMESS代码，笔者习惯将解压后得到的gamess/目录重命名，以便与未修改过的gamess区分开，例如

代码语言：javascript复制
mv gamess gamess_eda

若读者只使用一个版本，可以不重命名。进入目录，运行

代码语言：javascript复制
cd gamess_eda
../xeda-patch -v 2021-R2

这里笔者将xeda-patch脚本放在了上级目录，所以可以用../来运行。如果读者将其放在其他目录，则应写清楚路径。如果读者用的不是2021-R2，应在上述命令中替换相应的版本号。这种修改代码的方式与笔者的MOKIT程序里的modify_GMS1.sh脚本有相似之处。修改完之后，按正常编译GAMESS的步骤来即可，即先./config，根据屏幕提示填写信息，然后执行一系列make操作，细节可阅读之前的安装教程《GAMESS简易编译教程》。安装完后，笔者推荐在~/.bashrc中添加一条别名

代码语言：javascript复制
alias xeda=/home/$USER/software/gamess_xeda/rungms

这样我们就能在任意目录下直接使用xeda命令，而无需写一长串路径。

2. 示例

2.1 H2O…NH3复合物

我们在PBE0-D3BJ/aug-cc-pVDZ级别下做个GKS-EDA计算，示例如下

代码语言：javascript复制
 $CONTRL SCFTYP=RHF RUNTYP=EDA ICHARG=0 MULT=1 NOSYM=1 ISPHER=1
  DFTTYP=PBE0 $END
 $SYSTEM MWORDS=1000 $END
 $BASIS GBASIS=ACCD $END
 $DFT NRAD0=99 NLEB0=590 NRAD=99 NLEB=590 IDCVER=4 $END
 $LMOEDA EDATYP=GKS MATOM(1)=3,4 MCHARG(1)=0,0 MMULT(1)=1,1 $END
 $DATA
GKS-EDA of H2O...NH3 complex at PBE0-D3BJ/aug-cc-pVDZ
C1 1
 O   8.    0.16216217    0.62162161    0.00000000
 H   1.    1.12216217    0.62162161    0.00000000
 H   1.   -0.15829241    1.52655744    0.00000000
 N   7.    2.92216217    0.62162161    0.00000000
 H   1.    3.33225620   -0.28964503    0.03762964
 H   1.    3.13047339    1.11600642    0.84391356
 H   1.    3.29507043    1.12299101   -0.78074846
 $END

其中GBASIS=ACCD表示aug-cc-pVDZ基组，在计算量允许情况下，读者可采用cc-pVTZ或aug-cc-pVTZ基组获得更佳结果，此处仅为示例。$DFT中的格点设置是笔者的偏好，目的在于尽可能接近高斯DFT的ultrafine格点。现今可以不写其中的NRAD0和NLEB0，开发者告诉我现在程序会自动把NRAD和NLEB赋给NRAD0和NLEB0。IDCVER=4表示D3BJ。注意元素和坐标书写处有5列数据，比高斯gjf文件多出一列核电荷数（记得写小数点）。接着提交任务

代码语言：javascript复制
xeda h2o_nh3.inp 00 4 >& h2o_nh3.gms &

其中00表示GAMESS实际可执行程序gamess.00.x中的数字，4表示4核并行，输出文件为h2o_nh3.gms。如果你机器上默认的shell是/bin/dash而非常见的/usr/bin/bash，则应改为运行以下命令

代码语言：javascript复制
xeda h2o_nh3.inp 00 4 >h2o_nh3.gms 2>&1 &

注意无论任务正常结束与否，临时文件目录中会有残留文件，记得及时清理。在同一位置再次提交前，也应先清理临时文件。输出文件末尾几十行处有能量分解结果

代码语言：javascript复制
 -------------
 ALL BASIS SET                                       HARTREE          KCAL/MOL
 -------------
 ELECTROSTATIC ENERGY                  ES=         -0.024925            -15.64
 EXCHANGE ENERGY                       EX=         -0.040450            -25.38
 REPULSION ENERGY                     REP=          0.073507             46.13
 POLARIZATION ENERGY                  POL=         -0.010989             -6.90
 GRIMME DISP CORRECTION                DC=         -0.000830             -0.52
 ELECTRON CORRELATION                CORR=         -0.006362             -3.99
 TOTAL INTERACTION ENERGY               E=         -0.010049             -6.31

根据自己的需求，还可以在输入文件中添加隐式溶剂模型关键词。注意如果能量分解结果十分离谱，大概率是中间某一步计算的SCF没有收敛，GAMESS程序对于SCF不收敛是不会在末尾报错的（有点坑），只能自己翻看输出文件找报错，关于这点请往下看第3部分。

2.2 切断化学键得到的片段

上述示例是分子复合物，两个单体间没有键连关系。而GKS-EDA还支持分子内划分片段，不过要注意，每个片段的自旋要合理设置，例如将乙烷划分为两个甲基，两个基团都是二重态，在输入文件中要特地将两个片段的自旋多重度的数值写成相反数，以此保证整体Sz=0

代码语言：javascript复制
 $CONTRL SCFTYP=UHF RUNTYP=EDA ICHARG=0 MULT=1 NOSYM=1 ISPHER=1
  DFTTYP=PBE0 $END
 $SYSTEM MWORDS=1000 $END
 $BASIS GBASIS=ACCD $END
 $DFT NRAD0=99 NLEB0=590 NRAD=99 NLEB=590 $END
 $LMOEDA EDATYP=GKS MATOM(1)=4,4 MCHARG(1)=0,0 MMULT(1)=2,-2 $END
 $DATA
GKS-EDA of CH3-CH3 at UPBE0/aug-cc-pVDZ
C1 1
 C   6.   -1.72972983    2.31081078    0.00000000
 H   1.   -1.37307541    1.30200077    0.00000000
 H   1.   -1.37305699    2.81520897   -0.87365150
 H   1.   -2.79972983    2.31082396    0.00000000
 C   6.   -1.21638761    3.03676705    1.25740497
 H   1.   -0.14638763    3.03658448    1.25750243
 H   1.   -1.57288267    4.04563338    1.25730739
 H   1.   -1.57321987    2.53248124    2.13105625
 $END

由于这是化学键强度的结合能，所以数值会大很多

代码语言：javascript复制
 -------------
 ALL BASIS SET                                       HARTREE          KCAL/MOL
 -------------
 ELECTROSTATIC ENERGY                  ES=         -0.211366           -132.63
 EXCHANGE ENERGY                       EX=         -0.302565           -189.86
 REPULSION ENERGY                     REP=          0.631525            396.29
 POLARIZATION ENERGY                  POL=         -0.240778           -151.09
 ELECTRON CORRELATION                CORR=         -0.049101            -30.81
 TOTAL INTERACTION ENERGY               E=         -0.172285           -108.11

3. 用MOKIT自动产生GKS-EDA输入文件

实际研究复杂体系时，使用者经常碰到以下问题：

（1）GAMESS的SCF收敛性较差，收敛技巧较少，大量同学碰到SCF不收敛，不知如何解决或无法解决；而且GAMESS目前不支持检验UHF/UDFT波函数稳定性，对于复杂体系，即使通过调节各种关键词好不容易让SCF收敛了，波函数是否稳定也无法知道。

（2）GAMESS中的基组写法晦涩难懂，混合基组、赝势令人头疼。

（3）不同色散校正、隐式溶剂模型等关键词要逐个查阅GAMESS手册。

为了（暂时地）解决这些问题，笔者一年多前在MOKIT中添加了小程序frag_guess_wfn，已经有一些用户使用过了。它可以产生各种能量分解（SAPT0, MOROKUMA-EDA, GKS-EDA等）的输入文件，且含有必要的关键词、坐标、基组和轨道数据，这样在GAMESS里各个片段和整体的SCF都能极速收敛，且用户不需要懂GAMESS输入文件各种复杂格式，只需要会写高斯gjf文件即可。frag_guess_wfn小程序会自动调用高斯做HF/DFT计算，如果是UHF/UDFT计算，还会自动检验波函数稳定性，确保波函数稳定，然后将各片段的轨道写入inp文件。这里我们举一个[Ti(H2O)6]3 的例子，高斯gjf文件为

代码语言：javascript复制
%mem=56GB
%nprocshared=28
#p UPBE1PBE/gen guess(fragment=2)

{gks}

3 2 3 2 0 1
Ti(fragment=1)    0.00000000     0.00000000     0.00000000
 H(fragment=2)    0.00000000     2.68428109    -0.78751709 r
 H(fragment=2)    2.68128062     0.78751744     0.00000000
 H(fragment=2)    0.00000000    -2.68428109     0.78751709
 H(fragment=2)    0.78751709     0.00000000    -2.68228109
 H(fragment=2)   -2.68128062     0.78751744     0.00000000
 H(fragment=2)   -0.78751709     0.00000000     2.68228109
 O(fragment=2)    0.00000000     2.10100000     0.00000000
 H(fragment=2)    0.00000000     2.68428109     0.78751709
 O(fragment=2)    0.00000000     0.00000000     2.09900000
 H(fragment=2)    0.78751709     0.00000000     2.68228109
 O(fragment=2)    0.00000000     0.00000000    -2.09900000
 H(fragment=2)   -0.78751709     0.00000000    -2.68228109
 O(fragment=2)    2.09800000     0.00000000     0.00000000
 H(fragment=2)    2.68128062    -0.78751744     0.00000000
 O(fragment=2)   -2.09800000     0.00000000     0.00000000
 H(fragment=2)   -2.68128062    -0.78751744     0.00000000
 O(fragment=2)    0.00000000    -2.10100000     0.00000000
 H(fragment=2)    0.00000000    -2.68428109    -0.78751709

Ti 0
def2TZVP
****
H O 0
def2SVP
****

这实际上就是一种Gaussian做BSSE或片段组合波函数的文件格式。注意标题行要写{gks}，用来告诉frag_guess_wfn小程序产生GKS-EDA输入文件；第二个片段中的第一个原子坐标后写了一个字母r，用来告诉frag_guess_wfn对水分子做RPBE0计算，因为我们知道水分子做限制性闭壳层计算即可；而对Ti3 或复合物做UPBE0计算。这样既能节约计算时间，又不影响结果。有读者可能会问，那这个gjf文件里这么多fragment，有简便方法来书写么？笔者提供2种办法：（1）使用专业的文本编辑器（例如Sublime Text3，UltraEdit，VScode等）进行列操作；（2）使用GaussView主面板上的Tools -> Atom Groups定义片段，做过BSSE或片段组合波函数的读者对这个肯定不陌生。提交任务

代码语言：javascript复制
frag_guess_wfn Ti_H2O6.gjf >& Ti_H2O6.out &

大概7分钟后生成Ti_H2O6.inp文件，不用修改，直接提交给GAMESS即可

代码语言：javascript复制
xeda Ti_H2O6.inp 00 28 >& Ti_H2O6.gms &

其中所有的SCF步骤不超过4圈，很快就能算完。这个例子如果不借助frag_guess_wfn小程序的话、仅用GAMESS硬算的话，会有几个麻烦之处：

（1）需自己写GAMESS格式的混合基组；

（2）由于整体是开壳层计算，对水分子片段也只能做开壳层计算，但算出来其实是闭壳层结果，多花了机时；

（3）复合物SCF收敛困难，且波函数稳定性无法保证。

至于frag_guess_wfn小程序支持自动产生的其他类型能量分解文件，以后再介绍。

如何引用：

（1）使用XEDA模块做完GKS-EDA计算，发表正式文章时需引用相关文章，见

http://xacs.xmu.edu.cn/citation

（2）若使用了MOKIT中的frag_guess_wfn小程序自动产生输入文件，则应引用MOKIT，笔者提供了EndNote引用文件，见

https://gitlab.com/jxzou/mokit/-/tree/master/doc

若要以文本形式引用，见

https://gitlab.com/jxzou/mokit/-/blob/master/README.md

云函数 serverless 小程序 微信 游戏

0 人点赞