仓库源文，站点原文

layout: post title: GaussView使用cube文件作图时的成键问题 categories:

• 科 tags:
• 编程

2016-01-11 20:33:59

在使用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                3
NBond                                      I   N=          16
           1           2           1           1           2           1
           1           1           1           0           3           1
           1           2           1           1
IBond                                      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           0
RBond                                      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文件中成键信息内容. 这种方法不好的地方在于要重新做一个计算. 更好一些的办法当然是写一点代码来完成这个工作了. 丁迅雷的程序无法下载, 所以我不知道他的作法. 下面是我所用的脚本.

usage="\
>>>>>>>>>>>>>>>>     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=0
    while(Nblank<3) {
        getline txt
        gsub(/[[:space:]]/, "", txt)
        if(length(txt)==0) Nblank++
    }
    Natm=0
    getline
    while($1~/[1-9]/) {
        Natm=$1
        for(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=5
    print "MxBond                                     I                "MxBond

    print "NBond                                      I   N=          "Natm
    for(i=1; i<=int(Natm/6); i++) {
        for(k=1; k<=6; k++) {
            NBond=0
            for(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=0
        for(j=1; j<=Natm; j++) if(RBond[i,j]!=0) NBond++
        printf"%12d", NBond
    }
    print ""

    print "IBond                                      I   N=          "Natm*MxBond
    for(i=1; i<=Natm; i++) {
        NBond=0
        for(j=1; j<=Natm; j++) if(RBond[i,j]!=0) { NBond++; printf"%12d", j }
        for(j=1; j<=MxBond-NBond; j++) printf"%12d", 0
        print ""
    }

    print "RBond                                      R   N=          "Natm*MxBond
    for(i=1; i<=Natm; i++) {
        NBond=0
        for(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", 0
        print ""
    }
}
' $Fgjf >__Bond

awk '
    /Virial Ratio/ { print }
    /MxBond/, /Virial Ratio/ { while(getline < "__Bond") print; next }
    { print }
' $Fchk > ${Fchk%.*}~ReBond.fchk

rm -rf __Bond

这个脚本需要两个输入文件, 一个是含有成键信息的gjf文件, 一个是要修改成键信息的Fchk文件. 执行后, 脚本会将gjf文件中的成键信息写到Fchk文件中, 达到了我们的目的. 需要注意的是, 在这个脚本中, 我将MxBond设为了定值5以方便格式输出. 对一般的分子, 这个值足够了, 但如果有需要增大这个值, 相应的输出格式也要修改.