Gromacs的g_spatial产生的格点文件分析工具gmxgrid v1.2

下载地址/usr/uploads/file/20150605/20150605005449_28617.rar


gmxgrid V1.2
Convert the cube file generated by g_spatial to recognizable format at specific plane
Programmed by Sobereva, 2009-Oct-29

bug report: sobereva[at]sina.com

debut:2009-Feb-14
V1.1:2009-JUL-13,修复一个原点位置错误的bug
V1.2:2009-Oct-29,读取的变量双精度改为单精度,部分细节进行了改进,内存使用量减少为以前的1/4,增加了内存占用空间提示和读取进度提示。输出精度小数点后扩展为6位。

gromacs提供了g_spatial工具,可以生成gaussian格点文件来描述分子在各个位置的密度分布。虽然VMD也提供了volumeslice显示方式可以显示格点数据的切面图,gaussview4.1可以显示g_spatial生成的格点数据的自定义平面的等值线图,包括molden也可以实现同样功能。但是显示效果不好,毕竟是不是专业数据处理和绘图的软件,不适合用于文献的插图。
gmxgrid这个程序主要有三个功能,第一个是把格点文件的数据提取出来,输出每个数据点的坐标,即程序里的第1项。第二个功能是选取指定的XY/YZ/XZ平面,得到这个平面上的数据点以便做图,即程序里的第2-4项。第三个功能是将一定范围内平面的数据取平均输出,即程序里的第5-7项,比如Z值在30至33.5埃之间的平面的数据取平均输出。

gmxgrid的用处:例如要研究在分子平面中溶剂在分子周围分布情况,可以依如下方法进行:在VMD里旋转要研究的那个分子,使之与坐标轴组成的面平行或垂直,比如让分子内的六元环平行于X-Y平面,保存为新结构,用trjconv把轨迹向这个结构fit,用g_spatial得到溶剂密度的格点文件。用gmxgrid提取格点中那个分子所在平面的数据(或者z值在某个范围内的平面的密度和),用sigmaplot、origin做图,然后把分子图形套进那个图里(可以用相同方法做这个分子的密度分布,将两个图重叠,以确定分子在图中的位置)。


注意:格点文件中所用单位为波尔半径,等于0.529177249埃,读入VMD、gaussview后会被转换为埃。为了与可视化程序保持一致,本程序在读入格点文件时已经将单位转换为埃,输入参数和输出结果都以埃为单位。

读取特别大的格点文件,如几百兆,若程序刚载入文件就自动关闭,请确保有足够大的内存空间,应为总格点数*4*4字节。另外将虚拟内存设大。

关于g_spatial格点的单位,有两种情况,一种是默认情况,一种是加了-nodiv。默认情况格点文件里的数值都会除以一个归一化因子。加了-nodiv情况下,假设格点宽度就是-bin所设数值,即默认0.05nm,你计算了一个1ns的轨迹,得到其中一个格点的数值是1.8,就是说:在这1ns中,在这个0.05nm^3的空间内,平均有1.8个粒子。
粒子就是指的你选的group里的所有粒子。

以下关于高斯格点文件的介绍内容与gsgrid的readme.txt相同,但要注意g_spatial的格点文件与高斯生成的格点文件区别在于后者每6个数据换行一次,且最低维循环一遍后会换行,数据使用科学记数法。前者的每行内容就是最低维数据的一次循环,其间不换行,中维数据循环一次后会空一行。gsgrid与gmxgrid的区别在于gmxgrid处理的是g_spatial的文件格式,并添加了第5-7项功能。


格点文件介绍:
例如水的静电势的格点文件
Title Card Required potenial //顾名思义
Electrostatic potential from Total SCF Density //顾名思义
3 -4.970736 -4.970736 -4.761332 //原子数 原点的X/Y/Z坐标
80 0.125841 0.000000 0.000000 //第一个“坐标轴”上有80个数据点,每个数据点间隔为0.12584(以bohr为单位)
80 0.000000 0.125841 0.000000 //第二个“坐标轴”
80 0.000000 0.000000 0.125841 //第三个“坐标轴”
8 8.000000 0.000000 0.000000 0.209404 //原子序号(氧),电子数,X/Y/Z坐标(以bohr为单位)
1 1.000000 0.000000 1.481500 -0.837616
1 1.000000 0.000000 -1.481500 -0.837616
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=originx+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,彼此正交。