请选择 进入手机版 | 继续访问电脑版

分子模拟论坛 Molecular Simulation Forums

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 362|回复: 9

[原创]Gaussian格点文件分析、运算程序GsGrid V1.6

[复制链接]

83

主题

452

帖子

993

积分

高级会员

Rank: 4

积分
993
发表于 2008-9-25 23:26:00 | 显示全部楼层 |阅读模式
本帖最后由 sobereva 于 2010-3-16 18:28 编辑

最新版本请于官方网站下载

http://gsgrid.codeplex.com/




GsGrid V1.6
Extract data from Gaussian grid file and grid file calculation
Programmed by Sobereva, 2009-Oct-7
Bug report or recommend:
[color=]sobereva@sina.com
最新版本可来信索取或由我的blog提供的地址下载:
http://hi.baidu.com/sobereva/blog/item/605942a8dd21ffb6ca130c5b.html
1. 更新记录:
V1.0 最初版本,包含前四项功能。
V1.1 将输入输出单位由以前的波尔半径改为埃。加入了输出一定范围内平面的平均值的功能。
V1.3 修正了以前版本严重bug(即提取XZ平面数据得到的结果有时会偏差较大),新增功能8、功能9,用于得到自定义平面的格点数据。
V1.4 提高了output.txt中的数据输出精度(坐标精度不变),小数部分扩展为15位(整数部分最大5位),加入了格点文件计算功能。
V1.5 修正格点文件类型判断方式,使cube关键字生成的分子轨道格点文件可直接被gsgrid正确读取。新增支持多分子轨道格点文件。
V1.5.1 修正一个载入某些格点文件时出错的bug。
V1.5.2 修正读取多MO的格点文件失败的bug。
V1.6 加入了提取等值面上格点的功能,并可以将另外格点文件投影到某等值面上
2. GsGrid简介
Gaussian提供了cubegen工具,以及cube、prop选项可以得到格点文件。缺点是数据点没有对应的坐标信息,而且也无法提取某个平面的数据作图。虽然支持格点文件的可视化软件很多,GaussView较新版本及ChemCraft可以直接显示某个自定义平面上的等值线(前者支持)或填色图(后者支持),但是程序自建绘图功能极为有限,较为粗糙,不适合作为文献插图,无法导出数据在专业作图软件中自定义绘制,另外这两个软件也是收费的。至于molden,没有纯粹的windows版本,界面颇不友好。GsGrid主要为解决这些问题制作。
对于格点运算,虽然ChemCraft和cubman都能实现,但是所支持的运算操作各有不足,gsgrid从v1.4开始支持丰富的格点运算操作,范围超过ChemCraft和cubman支持操作的并集。
虽然很多软件,如果gaussview、Moliso等都可以将某个属性投影到某个等值面上并以颜色表示数值大小,比如最常见的就是将分子静电势(MEP)投影到电子密度为0.0001的等值面,用以研究分子VDW表面的静电势。然而这些软件并不提供将等值面上的数值输出的功能,只能用颜色显示,GsGrid弥补了这一不足。
2. GsGrid功能介绍
GsGrid目前支持6大类功能
第一类功能(程序中第1项):把格点文件的数据提取出来,附上每个数据点的坐标。见例1。
第二类功能(程序中第2-4项):选取指定的XY/YZ/XZ平面,得到这个平面上的数据点以便做图。见例2。
第三类功能(程序中第5-7项):将一定范围内平面的数据取平均输出,比如Z值在-3至3.5埃之间的平面的数据取平均输出。
第四类功能(程序中第8-9项):通过指定三个原子,或者自行输入三个点的坐标,定义一个平面,得到此平面的上的格点数据。此平面上的数据点可以直接输出,但是数据点并不在一个平面上,不方便做图。程序可以将此平面与XY平面的相交线做为轴,旋转此平面,使数据转移至XY平面,即此时数据点Z值皆为0,随后可通过绘图软件做contour图。见例3。
第五类功能(程序中第10项):
格点文件中每个数据对某个常数进行加减乘除(子菜单1/3/5/7项)
格点文件中每个数据对另一个格点文件对应项进行加减乘除(子菜单2/4/6/8项)
格点文件中每个数据求幂(子菜单第9项);格点文件中每个数据对另一个格点文件对应项进行平方和、平方差、求平均操作(子菜单10/11/12项)
格点文件中每个数据求绝对值(子菜单第13项)
操作完毕后会输出新的格点文件,新的格点文件可继续使用gsgrid进行操作。见例4、例5。
第六类功能(程序中第11-12项):第11项功能为提取某个等值面上的格点。第12项功能能够将代表某个属性的格点文件A投影到格点文件B的数值所定义的等值面上,输出文件中包含这些等值面上的格点的坐标和A格点文件中坐标相对应的格点的数值。见例6、例7。
3. 例子
以下例子建议初学者从头看一遍。
例1:得到空间中每个坐标点的静电势
用高斯计算某个体系,写明%chk,然后用formchk将chk转换为fch,例如得到marioholic.fch。cubegen.exe在高斯目录下,将cubegen程序与marioholic.fch放到同一个文件夹,在dos中进入相应目录下运行:
cubegen 0 potential marioholic.fch aneta_langerova.cub 0 h
就得到了静电势格点文件aneta_langerova.cub。
(也可以在这个文件夹下用文本编辑器建立一个批处理文件,比如a.bat,内容就是上面这条命令,然后双击运行a.bat)
运行gsgrid,比如输入:(注:下文//后面皆代表注释,不要输入。换行处就是敲回车)
c:\aneta_langerova.cub //若.cub文件与gsgrid在同目录,只需输入文件名)
1 //功能1
即在相同目录下得到output.txt。里面包含aneta_langerova.cub里每个数据点的数值及其坐标。
------------------------------------
例2:画Z=0平面的电子密度图
先得到电子密度格点文件,运行:cubegen 0 density mizuki_nana.fch otoboku.cub -4 h
运行gsgrid,输入:
c:\otoboku.cub
2 //功能2,提取某XY平面的数据
0 //设定此XY平面的Z值为0
得到output.txt,和例1一样。由于格点文件数据是离散的,故这里面的数据点只包含.cub中的Z值最接近于0的那个平面中的数据点。然后可以放到sigmaplot等软件做图。以苯分子为例子(压缩包内c6h6.fch),以此法作图得到benzene-density.PNG及benzene-density-contour.GIF
再比如想得到X在-3.3至5.2范围内的YZ平面上的数据点的平均值,可使用程序第6个功能,运行gsgrid,输入
c:\otoboku.cub
6 //功能6
-3.3 5.2 //Z值范围
------------------------------------
例3:得到三个原子定义的平面的静电势图
设某个分子有三个原子编号为3、6、7,想得到这三个原子定义的平面上的静电势图,但是这三个原子并不平行于XY/YZ/XZ,故不能用程序的2、3、4项功能通过例2的方法做图。在V1.3版中添加了任意平面做图可以实现这个目的。
cubegen 0 potential ethane.fch 1.cub -4 h //ethane.fch在附件中,这里使用fine格点,否则图像会破碎,见后面说明。
运行gsgrid,输入:
c:\1.cub
8 //功能8
3,6,7 //输入三个原子序号
0.014 //以0.014埃为格点-平面临界距离
y //将平面数据通过旋转投影到XY平面上
这样得到的结果就像例2,数据点z值皆为0,可以直接用绘图软件做图,使用sigmaplot得到的结果如附件中ethane-PES.PNG所示。
再比如想自定义三个点得到定义的平面的数据,运行gsgrid,输入:
c:\1.cub
9
2.3,1.1,2.55 //第1个点的XYZ
0.23,1.33,4.5 //第2个点的XYZ
0,-1.29,3.11 //第3个点的XYZ
0 //用默认的格点-平面临界距离
n //不投影至XY平面
------------------------------------
例4:得到某体系第4个分子轨道在某XY平面的密度分布
首先需要得到分子轨道波函数格点文件,运行cubegen 0 mo=4 divokej_bill.fch 1.cub 0 h。
求每个格点数据的2次方,因为波函数平方为概率密度。运行gsgrid,输入:
c:\1.cub
10 //第10项功能,即计算功能
9 //子菜单第8项功能,求幂
2 //设定指数,输入的数除了整数外,也可以是小数、负数。要注意底数是负数时指数须为整数,否则结果为NaN。
z.cub //输出的格点文件,不写路径则生成在当前路径下
假设这个分子轨道由两个电子占据,故每个格点数据应再乘以2,再次启动gsgrid,输入:
z.cub
10 //第10项功能,即计算功能
5 //子菜单第5项功能,乘以一个常数
2 //每个格点数据乘以2
lusaint.cub //输出的格点文件
最后再启动gsgrid,提取lusaint.cub的某XY平面数据即可,见例1。
再例如我们想得到mo=4和mo=5两个分子轨道上占据的电子密度和的格点文件。mo=4的格点文件仍是上面的1.cub,mo=5的格点文件是2.cub。运行gsgrid
1.cub //先输入第一个格点文件
10 //第10项功能,即计算功能
10 //子菜单第10项功能,求平方和,即A^2+B^2
2.cub //第二个格点文件
Bakuretsu_Tenishi.cub //输出的格点文件
若每个轨道是双占据的,再把Bakuretsu_Tenishi.cub用gsgrid乘以常数2即可。
实际上用gsgrid把1.cub和2.cub单独求平方,然后再相加得到的结果和上面一样,这里用子菜单第10项功能是因为比较方便,只需要一步而不必三步。
要注意,两个格点文件间运算,两个格点文件必须是相同体系,有相同的格点数目。
------------------------------------
例5:多分子轨道格点文件的读取
每个高斯格点文件可以包含多个轨道,从gsgrid V1.5开始支持这样的格点文件。这里仍然使用例4中的第二个例子,需要获得MO=4和MO=5的总电子密度格点文件。为了方便,我们现在把mo=4和mo=5的格点都写进一个文件里,例如cubegen 0 mo=4,5 bitboys.fch sob.cub 0 h。运行gsgrid,输入
sob.cub
1 //选取的格点文件若包含多个轨道信息,则会出现选择轨道的提示。这里输入1代表选择第一条轨道,即MO=4
10 //第10项功能,即计算功能
10 //平方和运算
sob.cub //第二个格点文件,因为MO=5的轨道就在sob.cub中,所以仍用sob.cub
2 //选格点文件中第二条轨道,即MO=5
z.cub //输出文件
闭壳层体系时再把z.cub的内容乘2即可。这样将多条轨道合并在一个格点文件中用起来比较方便,但是会造成格点文件成倍增大,gsgrid读取也会变慢。
第一个格点文件和第二个格点文件可以包含的轨道数目不同,也可以不是一类内容,如分子轨道与电子密度值,只要格点在数目、空间位置上是对应的即可。
使用gsgrid也可以提取多分子轨道格点文件中某一个轨道的格点文件,即选择那条轨道,使用比如加0或者乘以1的运算。
高斯的分子轨道(包括多分子轨道)的格点文件与其它类型格点文件的主要区别在于前者的原子数为负,在原子坐标后面会加入分子轨道信息(见后面附带的格点文件格式说明)。如果第一个格点文件是分子轨道格点文件,则运算后生成的新格点文件也将是分子轨道类型。若格点文件中原子数为正/负,则必须无/有分子轨道信息,否则格点文件无法被gsgrid正确处理。
------------------------------------
例6:提取某个等值面的数据点
这个例子要用到1.6加入的功能11。下面我们要得到自带的test文件夹下c6h6.fch体系(苯,3-21G)的第19个分子轨道波函数的绝对值为0.02的等值面的格点。首先得到第19个分子轨道的格点文件:
cubegen 0 mo=19 c6h6.fch c6h6-mo-19.cub -3 h
使用功能11、12时最好使用中等及以上精细度的格点文件,所以这里使用-3参数即中等精细度。
运行gsgrid,依次输入:
c6h6-mo-19.cub
11//第11项功能,提取某个等值面的格点
0.02,0.02//这两个数分别定义下界和上届的数值,如果某格点的数值在这个范围内,则这个格点就会输出到output.txt。这里上界和下届输入的数值一样,则程序会将这个数值的±3%作为上下界,即提取0.0194至0.0206范围内的格点(若不允许有±3%的偏差,得到的数据点太少)。
现在得到了output.txt,但这里只有分子轨道波函数为正值的部分,我们还想得到等值面为-0.02的格点,故把刚才得到的output.txt备份,再次在gsgrid中运行上述命令,0.02,0.02替换为-0.02,-0.02即可。
现在我们将这些格点放到origin里面做三维散点图,正值与负值部分分别用红色和绿色表示,得到了isosurface.PNG文件的左图。可以看到,我们提取的等值面上的格点是正确的,形状与gaussview中显示的此分子轨道波函数的0.02的等值面完全一致。
------------------------------------
例7:将静电势投影到电子密度为0.001的等值面
首先分别得到电子密度和静电势的格点文件,注意它们必须拥有相同格点数,每个格点坐标相等。
cubegen 0 density c6h6.fch c6h6-m-density.cub -3 h
cubegen 0 density c6h6.fch c6h6-m-eps.cub -3 h
运行gsgrid,依次输入:
c6h6-m-density.cub//通过这个格点文件确定等值面的格点
12//选功能12
0.001//等值面数值为0.001
4//设定上下限偏差,用百分数表示。这里即设定数值在0.001±4%范围内的格点看做等值面格点
c6h6-m-eps.cub//将这个静电势格点文件投影到刚才设的等值面上并输出
这样我们得到的output.txt中的数据点的坐标就是c6h6-m-density.cub文件中数值为0.001±4%范围所定义的等值面上的格点的坐标,而这些点的数值就是c6h6-m-eps.cub文件中坐标相对应的点的数值。
------------------------------------
关于功能8和功能9的细节说明:
由于格点文件的数据点在空间上是规律排布的,若自行定义一个平面,数据点显然不会恰好在平面上。所以想得到某个平面上的数据,就必须将平面附近的数据点投影到平面上。在gsgrid中,只有数据点与平面垂直距离小于“格点-平面临界距离”的数据点才会被投影,符合要求的数据点投影的方法是沿着法线方向移动至平面。格点-平面临界距离的设定会影响实际效果,一般用默认值即可(输入0使用建议值)。
有时通过这种投影方法获得的某平面的数据在做等值线图时会发现图像虽然形状正常,但略有破碎,或者个别点数据怪异。此时可降低格点-平面临界距离来解决,例如改为默认值的1/2、1/4再重新做图,但是格点-平面临界距离如果太近的话,投影到平面上的格点数太少,图像会变得粗糙。最好的解决方法是生成格点文件时使用fine格点,此时格点密度加大,投影到平面上的数据准确度较好,可消除图像破碎问题,此时建议将格点-平面临界距离设为推荐值的1/2。
有时得到的等值线图有棱有角,是因为格点文件不够精细造成,亦可通过使用fine格点解决,也可以通过改变等值线或色彩范围、调节格点-平面临界距离来得到一定改善。
将平面的数据转移至XY平面的过程,是以自定义平面与XY平面的交线为轴,旋转自定义平面来实现的。不会使平面上数据得到的图像有任何变形,相当于眼睛垂直于自定义平面观看。
如果定义的平面与XY/YZ/XZ平面平行,程序会予以提示,此时请分别用2、3、4项功能代替,以得到更好结果。
关于8、9项功能使用若有问题或结果异常,请E-mail来信。
4. 使用注意
使用cubegen时,末尾一定要用h参数。
cubegen倒数第二个参数可以是-2,-3 和-4 分别对应于关键字Coarse,Medium和Fine。用了fine文件会很大,gsgrid读取也比较慢,但相应地数据点更精细,建议绘制文献图时使用fine。若设为0用默认值,大约与Medium相当。
gsgrid目前只支持X、Y、Z平移向量与坐标轴X、Y、Z平行的格点文件,即一般情况下的格点文件。
格点文件中所用单位为波尔半径,等于0.529177249埃,读入VMD、gaussview后会被转换为埃。为了与可视化程序保持一致,本程序在读入格点文件时已经将单位转换为埃,输入参数和输出结果都以埃为单位。
5. 使用技巧:
若想使程序以silent模式运行,也就是敲一个命令即完成指定任务,不必每一步都手动输入,可以编辑一个文件,比如batch.txt,内容如下:
1.cub
8
3,6,7
0.014
y
[空行]
[空行]
这样只要运行gsgrid < batch.txt,就可以完成例3中的任务,每一步都根据batch.txt的内容自动输入了。如果处理很多格点文件,这样会方便很多。
6. 附录
Gaussian格点文件格式简介:
例如水的静电势的格点文件
Title Card Required potenial//Title,不被gsgrid处理
Electrostatic potential from Total SCF Density //Title,不被gsgrid处理
3-4.970736-4.970736-4.761332//原子数(如果是分子轨道格点文件原子数为负值)原点的X/Y/Z坐标
800.1258410.0000000.000000//第一个“坐标轴”上有80个数据点,每个数据点间隔为0.12584波尔半径
800.0000000.1258410.000000//第二个“坐标轴”
800.0000000.0000000.125841//第三个“坐标轴”
88.0000000.0000000.0000000.209404//原子序号(氧),电子数,X/Y/Z坐标
11.0000000.0000001.481500-0.837616
11.0000000.000000-1.481500-0.837616
//若是分子轨道格点文件,这里有分子轨道信息。第一个数字为格点文件包含的分子轨道的数目,接下来是每个分子轨道的序号。例如3 1 5 7就是代表格点文件含3个分子轨道,分别是1、5、7号。分子轨道信息内容中每10个数字换一行。
7.33384E-03 7.33602E-03 7.33151E-03 7.31988E-03 7.30070E-03 7.27354E-03//每个点的数据,每六个换一行,E13.5格式。如果数据点数不是6的倍数,最低维循环后不足6个数据位置的地方会留空。
7.23796E-03 7.19353E-03 7.13981E-03 7.07636E-03 7.00279E-03 6.91867E-03
6.82364E-03 6.71732E-03 6.59939E-03 6.46955E-03 6.32753E-03 6.17312E-03
......
这些数据点数据输出的循环次序是:最低维=第三个“坐标轴”(此例即是Z轴),中维=第二个“坐标轴”,最高维=第一个“坐标轴”。
数据点与其对应的坐标关系是:
a(i,j,k)%x=originx+trans1x*(i-1)+trans2x*(j-1)+trans3x*(k-1)
a(i,j,k)%y=originy+trans1y*(i-1)+trans2y*(j-1)+trans3y*(k-1)
a(i,j,k)%z=originz+trans1z*(i-1)+trans2z*(j-1)+trans3z*(k-1)
trans1/x/y/z代表第一个“坐标轴”三个分量,(或曰:格点的平移矢量)
trans2/x/y/z代表第二个“坐标轴”三个分量
trans3/x/y/z代表第三个“坐标轴”三个分量
originx/y/z代表原点位置
i、j、k代表某点在各个“坐标轴”上的数据点序号。
如果用cubegen默认设置,一、二、三号坐标轴实际上就相当于笛卡尔坐标轴X/Y/Z,彼此正交。
如果手动设置网格,可以自定义“坐标轴”的方向,不平行于笛卡尔坐标轴,但这种情况gsgrid只能提取全部数据点(第一个功能),不能提取指定平面的数据点。
也可以用cube、prop关键字生成格点文件,参见高斯手册,格式与cubegen生成的一样,但更建议使用cubegen来做。




回复

使用道具 举报

10

主题

63

帖子

140

积分

注册会员

Rank: 2

积分
140
发表于 2008-9-29 15:09:00 | 显示全部楼层
sobereva,你太强了,thanks for the share

回复

使用道具 举报

10

主题

63

帖子

140

积分

注册会员

Rank: 2

积分
140
发表于 2008-9-29 15:15:00 | 显示全部楼层
ms用gaussian 本身也可以,然后用gaussian view 看,easy for beigner,也可以用moden来看 MO density
回复

使用道具 举报

83

主题

452

帖子

993

积分

高级会员

Rank: 4

积分
993
 楼主| 发表于 2008-9-30 04:13:00 | 显示全部楼层
gsview等不少软件也支持cub,不过显示的往往是三维表面,没法看二维截面图,需要做二维图不方便。
molden可以,只是我不太喜欢molden。所以写这个小程序,导出平面数据点用专业数据处理软件做图比较舒服,这样等值线刻度、样式等等都可以方便地自行调整。
回复

使用道具 举报

8

主题

48

帖子

110

积分

注册会员

Rank: 2

积分
110
发表于 2008-10-2 11:19:00 | 显示全部楼层
哈哈,确实确实
同意楼上说法
哈哈,我经常计算大的体系的时候很难得到chk文件,或者chk文件太大,不方便保存,我自己写了个程序,直接从out文件当中计算cube文件。或许我也可以这样干一把
回复

使用道具 举报

17

主题

67

帖子

155

积分

注册会员

Rank: 2

积分
155
发表于 2008-10-6 20:08:00 | 显示全部楼层
不太懂, 觉得东西不错, 下载了. 顶一个!
回复

使用道具 举报

18

主题

124

帖子

272

积分

中级会员

Rank: 3Rank: 3

积分
272
发表于 2008-10-17 02:43:00 | 显示全部楼层
嗯, 下载学习. 辛苦了!
回复

使用道具 举报

50

主题

906

帖子

1868

积分

金牌会员

Rank: 6Rank: 6

积分
1868
发表于 2009-2-3 20:27:00 | 显示全部楼层
强啊~~
回复

使用道具 举报

131

主题

447

帖子

1031

积分

金牌会员

Rank: 6Rank: 6

积分
1031
发表于 2009-2-21 10:31:00 | 显示全部楼层
我这里有一个.cube文件,可是当我运行才c:\*.cube 时,生成的output是0K大小。里面什么也没有呀。这是怎么一回事。另外我的*.cube文件不是gaussian生成的,用这个小程序行吗?
回复

使用道具 举报

1

主题

7

帖子

17

积分

新手上路

Rank: 1

积分
17
发表于 2009-2-27 11:07:00 | 显示全部楼层
下载了,谢谢楼主分享。
回复

使用道具 举报

Archiver|手机版|小黑屋|分子模拟论坛  

GMT+8, 2019-9-18 19:42 , Processed in 0.160866 second(s), 30 queries .

Powered by Discuz! X3.3

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表