博文
GaussView使用cube文件作图时的成键问题
|
在使用GaussView打开cube文件做等值面之类的图时, GaussView会自动根据原子坐标判断原子之间的连接性, 并以成键的方式显示出来. 由于cube文件中没有储存原子间的连接信息, 所以判断是根据GaussView设定的原子半径进行的. 可惜的是在GaussView中这些原子半径都是固定的, 我们没有办法改变. 更讨厌的是, 打开cube文件后, 如果你选中两个原子并改变了它们之间的成键情况, 这时Results菜单下的Surfaces/Contours功能就无法使用了. 这实在是非常不方便, 因为有些情况下GaussView给出的成键情况不能满足我们的需要, 在使用cube文件做图的时候我们需要改变成键情况.
网上的解决办法大致有两种, 一种是换用其他软件, 如VMD, Chemcraft来作图, 如这里的建议. 这不是我们要讨论的方法. 另一种方法则是丁迅雷给出的方法, 并给出了一个修改Fchk文件中成键信息的程序, 但语焉不详, 程序也无法下载. 为此, 我对第二种方法略加说明, 供大家参考.
解决这个问题的办法是同时使用Fchk文件和cube文件, 因为除cube文件外, 只有在打开Fchk文件时, Results -> Surfaces/Contours功能才可使用, 且还能导入cube文件作图, 而Fchk文件中是可以保存成键信息的. 因此, 我们只要先在Fchk文件中设定好了成键情况, 使用GaussView打开它, 再导入cube文件作图就能达到目的了.
所以, 现在的问题就是看看Fchk中是如何保存原子成键信息的. 使用GaussView打开一个测试文件, 成键情况如下:
查看Fchk文件中的内容, 发现成键信息保存在下面的文本中
MxBond I 3NBond I N= 16 1 2 1 1 2 1 1 1 1 0 3 1 1 2 1 1IBond I N= 48 11 0 0 3 4 0 2 0 0 2 0 0 6 7 0 5 0 0 5 0 0 9 0 0 8 0 0 0 0 0 1 12 13 11 0 0 11 0 0 15 16 0 14 0 0 14 0 0RBond R N= 48 1.00000000E+00 0.00000000E+00 0.00000000E+00 1.00000000E+00 1.00000000E+00 0.00000000E+00 1.00000000E+00 0.00000000E+00 0.00000000E+00 1.00000000E+00 0.00000000E+00 0.00000000E+00 1.00000000E+00 1.00000000E+00 0.00000000E+00 1.00000000E+00 0.00000000E+00 0.00000000E+00 1.00000000E+00 0.00000000E+00 0.00000000E+00 1.00000000E+00 0.00000000E+00 0.00000000E+00 1.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 1.00000000E+00 0.00000000E+00 0.00000000E+00 1.00000000E+00 0.00000000E+00 0.00000000E+00 1.00000000E+00 1.00000000E+00 0.00000000E+00 1.00000000E+00 0.00000000E+00 0.00000000E+00 1.00000000E+00 0.00000000E+00 0.00000000E+00对照一下容易发现:
MxBond为所有原子的最大连接数.
对上面的体系11号原子与1, 12, 13号原子都相连, 具有最大的连接数3.NBond相当于一个数组, 维数为原子个数, 保存了每个原子的连接数.
第一个值为1说明只有一个原子与1号原子相连.IBond相当于一个二维数组, 维数为原子个数×MxBond, 保存了与每个原子相连的其他原子的编号.
上面的体系共16个原子, Mxbond为3, 所以IBond的元素数目为16×3=48.
前三个值11 0 0说明和1号原子相连的三个原子的编号分别为11, 0, 0, 其中0表示空, 所以实际只有1个原子与1号原子相连, 这和NBond的数据一致.RBond是与IBond类似的二维数组, 但保存的是原子之间的键级, 用于显示. 0为不成键, 1为单键, 2为双键, 3为三键, 也可以取0.5, 1.5等数值.
第一值是1, 所以说明和1号原子所成的第一根键是单键, 这个图上的显示一致.
理解了这种格式, 我们就可以自己来修改这些信息了. 手动修改还是比较麻烦的, 只适用于比较简单的情况. 另外一种办法是先使用GaussView将成键情况修改好, 另存为gjf文件, 并保留geom=connectivity关键词和成键信息, 使用最低水平的方法做一个单点能计算, 获得Fchk文件, 然后使用这个Fchk文件中的成键信息替换掉原Fchk文件中成键信息内容. 这种方法不好的地方在于要重新做一个计算. 更好一些的办法当然是写一点代码来完成这个工作了. 丁迅雷的程序无法下载, 所以我不知道他的作法. 下面是我所用的脚本.
# Language: bashusage=">>>>>>>>>>>>>>>> ReBondFchk <<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>> Jerkwin@gmail.com <<<<<<<<<<<<<<<<>>>>>>>>>> 2016-01-06 20:31:08 <<<<<<<<<>> Usage: `basename $0` File.gjf File.fchk"[[ $# -lt 1 ]] && { echo "$usage"; exit; }[[$1=~.gjf ]]&&{ Fgjf=$1; Fchk=$2;}[[!$1=~.gjf ]]&&{ Fgjf=$2; Fchk=$1;}awk '/#/ { Nblank=0while(Nblank<3){getline txtgsub(/[[:space:]]/,"", txt)if(length(txt)==0) Nblank++}Natm=0getlinewhile($1~/[1-9]/){Natm=$1for(i=2; i<=NF; i+=2) RBond[Natm,$i]=$(i+1)getline}for(i=1; i<=Natm; i++){for(j=1; j<=i; j++){RBond[i,j]=RBond[j,i]}}MxBond=5print "MxBond I "MxBondprint "NBond I N= "Natmfor(i=1; i<=int(Natm/6); i++){for(k=1; k<=6; k++){NBond=0for(j=1; j<=Natm; j++)if(RBond[6*(i-1)+k,j]!=0) NBond++printf"%12d", NBond}print ""}for(i=6*int(Natm/6)+1; i<=Natm; i++){NBond=0for(j=1; j<=Natm; j++)if(RBond[i,j]!=0) NBond++printf"%12d", NBond}print ""print "IBond I N= "Natm*MxBondfor(i=1; i<=Natm; i++){NBond=0for(j=1; j<=Natm; j++)if(RBond[i,j]!=0){ NBond++; printf"%12d", j }for(j=1; j<=MxBond-NBond; j++) printf"%12d",0print ""}print "RBond R N= "Natm*MxBondfor(i=1; i<=Natm; i++){NBond=0for(j=1; j<=Natm; j++)if(RBond[i,j]!=0){ NBond++; printf"%16.8E", RBond[i,j]}for(j=1; j<=MxBond-NBond; j++) printf"%16.8E",0print ""}}' $Fgjf >__Bondawk '/Virial Ratio/{ print }/MxBond/,/Virial Ratio/{while(getline <"__Bond") print; next }{ print }' $Fchk>${Fchk%.*}~ReBond.fchkrm -rf __Bond这个脚本需要两个输入文件, 一个是含有成键信息的gjf文件, 一个是要修改成键信息的Fchk文件. 执行后, 脚本会将gjf文件中的成键信息写到Fchk文件中, 达到了我们的目的. 需要注意的是, 在这个脚本中, 我将MxBond设为了定值5以方便格式输出. 对一般的分子, 这个值足够了, 但如果有需要增大这个值, 相应的输出格式也要修改.
◆本文地址: http://jerkwin.github.io/2016/01/11/GaussView使用cube文件作图时的成键问题/, 转载请注明◆
https://blog.sciencenet.cn/blog-548663-949826.html
上一篇：单层螺线管周围的磁场
下一篇：构建表面修饰长链分子的纳米颗粒模型
