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，彼此正交。