谈谈体数据2:体数据格式转换工具vodaconv的使用及原理

谈谈体数据2:体数据格式转换工具vodaconv的使用及原理

文/Sobereva @北京科音
First release: 2012-Feb-14  Last update: 2012-Feb-22


前言:在《谈谈体数据1:介绍体数据》(http://sobereva.com/127)一文中寡人介绍了体数据的相关概念。raw是广为通用的体数据储存格式,可惜计算化学领域的可视化软件极少有支持raw格式的;而通用的渲染体数据的程序也没有能支持计算化学领域通用的Gaussian型cube文件(.cub)的。为了能够让体数据互用,架起领域间的桥梁,寡人开发了vodaconv程序用于格式转换,将在本文介绍。在《谈谈体数据》系列的接下来的文章中将会用到这个程序对raw和cube型体数据进行相互转换。同时,本文也将仔细介绍此程序的原理,将有助于读者了解与体数据相关的一些细节问题。

1 vodaconv程序介绍

体数据转换工具vodaconv支持的输入格式包括raw、Gaussian型cube文件(.cub)、OpenDX的格点数据格式(.dx)、Dmol3的.grd格式。支持的输出格式包括raw和cube。在未来的版本中将会支持更多格式。此程序只支持立方格子的格点文件,即平移向量平行于笛卡尔坐标轴。

此程序可以在此下载:/usr/uploads/file/20150611/20150611014131_45262.rar

程序会根据扩展名判断格式。如果输入的是.raw文件,需要由用户自行输入各个方向的格点数和格点间隔、起始位置以及数据类型。如果输入不对的话将导致此程序在读取.raw文件中自动退出。可以事先自行算一下输入的格点数和数据类型与当前文件是否相符,比如含有256*256*128个格点的16bit的raw文件,文件大小应当为256*256*128*16/8=16777216字节,如果文件大小与之不符,说明数据类型或者格点数与此文件不符。

程序目前支持的raw数据类型包括8bit整数、16bit整数和单精度浮点数三种。读入、输出raw文件时整数一律按照无符号整数形式处理,即8bit整数数据范围为0~255,16bit整数数据范围为0~65535。vodaconv.exe读写的raw文件是按照Big Endian顺序,而vodaconv_LE.exe是按照Little Endian顺序(含义将在后文解释)。

在读入格点数据后会提示是否对数据进行映射,需要先输入原先的数据范围上下限,超过上限的数据将被设为上限值,超过下限的数据将被设为下限值。然后需要输入映射后的数据范围上下限。原先的数据将按照上下限范围的比例映射到新的数值范围内。这个功能主要用于调节数值范围,以免以8/16bit整数方式输出raw文件时数据越界。利用这个功能也可以将16bit整数的raw格式转化为8bit整数的raw格式,即downsample过程。一个应用例子是,比如想把电子定域化函数(ELF)的cube文件转化成16bit整数raw文件,由于已知ELF函数值域是[0,1],原先数据范围就可以输入0,1,而映射后的数据范围就可以输入0,65535,这样ELF函数值域区间就均匀展开到了无符号16bit整数的记录范围里了。

之后会提示用户输入要输出的文件名,也是根据后缀名判断输出格式。在输出cube文件时,若之前载入的是.raw/.dx/.grd格式,由于其中没有定义原子信息,为了避免一些程序在读取cube文件时必须要求至少有一个原子造成的不兼容,此程序将自动在cube文件内添加一个位于原点的氢。


2 vodaconv程序的编译方法

程序压缩包内已包含win32平台下由Intel visual fortran 12.0预编译好的可执行文件。同时附带了源代码文件。

若想自行编译源代码,在Linux下执行ifort vodaconv.f90 -o vodaconv即可。不要妄图用gfortran编译,因为此程序用了integer*1的数据类型,gfortran不支持。

如果想在Windows下编译,用Compaq visual fortran或Intel visual fortran都可以。首先将vodaconv.f90拖进新建的项目里,然后设置stack size为较大值,比如256MB(否则由于默认栈内存过小,在载入较大的格点文件时会导致栈内存溢出),之后照常编译即可。在CVF中设置stack size的方法是选择"Project"-"Settings..."-"Link",在"Category"里选择"Output",在"Reserve"里输入"0x10000000"(即256MB,这也是CVF支持的上限)。在IVF中设置stack size的方法是选择“项目”-“XX属性”,选"Linker",将"System"里的"Stack Reserve Size"设为"268435456"(即256MB。当然也可以设得更大使程序能载入巨型格点数据)。

如果你希望你所编译出的本程序使用Little Endian顺序读写二进制文件,那么就按上面的步骤编译就行了,编译器默认就是用Little Endian。如果希望使用Big Endian顺序,那么在Linux下编译时应该加上-convert big_endian选项。如果是在IVF中,应当在项目的属性的Fortran-Compatibility中将Unformatted file conversion设成Big Endian。如果是在CVF中,进入"Project"-"Settings-"Fortran",选Compatibility分类,将Unformatted file conversion设成Big Endian。


3 关于含符号整数和无符号整数

这里介绍一下含符号和无符号整数与raw文件的关系。这里讨论的都是Fortran语言涉及的整数。

计算机程序中用到的整数数据一般是双字节整数(16bit整数。也叫短整数,short integer)、四字节整数(32bit整数。也叫长整数,long integer),在现代的编译器中(比如版本较新的ifort)还支持八字节整数。单字节整数,即8bit整数一般很少用到,因为它的数值范围太窄,一般没什么用。而且它并不是Fortran标准中定义的数据类型,但CVF、ifort等编译器是支持它的(用integer*1定义。在gfortran中不支持)。然而8bit整数却被普遍被raw文件用于记录低精度数据,因此我们不能忽视它。

无符号整数意思是整数值只为正,因此8bit无符号整数的范围是0~2^8-1,即0~255;16bit无符号整数是0~2^16-1,即0~65535。含符号整数则将一半记录空间用于记录负值,8bit和16bit含符号整数的数值范围分别是-128~127和-32768~32767。

在Fortran标准中只定义了含符号整数,然而raw格式中的整数数据都是使用无符号整数来记录的,所以必须搞清楚含符号和无符号整数间的转换关系,才能用Fortran程序正确载入整数型raw格式数据。

在计算机的内存或二进制文件中,四个十六进制数由小到大对应的16bit整数数值如下所示,括号中第一个对应的是含符号型整数值,第二个是无符号型整数值:
0x0000(0,0)->0x0001(1,1)->...->0x7FFF(32767,32767)->0x8000(-32768,32768)->0x8001(-32767,32769)->...->0xFFFF(-1,65535)。
即曰:对于0~32767的部分,含符号和无符号的整数是以相同方式记录的。对于含符号整型,如果大于了32767就会成为负值,此时这个负值加上65536就可以变为相应的无符号整型的数值。这个规则对于其它长度的整型数据也是一样。

由于Fortran用的都是含符号整型,这也就是说,raw文件记录的8bit/16bit无符号整型数据会被当作是含符号整型数据被读入,因此需要在读入后把所有负值数据都加上256(对于8bit)和65536(对于16bit),这样才是正确的数值。

在通过write语句输出8bit和16bit整数型raw文件前,需要先分别利用int1()和int2()函数将内存中的数据转换为8bit和16bit含符号整型。对于8bit整型,大于127的数据将被自动减去256成为负数使之能被含符号整型变量表示;对于16bit则会自动减去65536。这个过程和载入时需要做的变换是互逆的。


4 二进制文件的字节记录顺序问题

不同的平台、不同编译器设定下产生的不同程序所利用、生成的二进制文件中的多字节数据有不同的字节记录顺序。8bit整数就一个字节所以无需考虑此问题,而16bit整数有两个字节,就会牵扯到这个问题。简要来说,例如32444这个16bit整数对应的十六进制数是0x7EBC,这个数据在二进制文件中如果低字节记录的是0xBC,高字节是0x7E,那么就称为Little Endian(小端)顺序,用Ultraedit等支持十六进制的编辑器打开此文件就会看到相应位置显示的是BC 7E。若低、高字节分别记录的是0x7E和0xBC,就称为Big Endian(大端)顺序。Windows、Linux等系统下的程序通常用的是Little Endian,而据说Mac OS用的是Big Endian。

如果你发现用vodaconv.exe(用的是Big Endian)转换的16bit整数型raw文件在载入其它一些体数据可视化程序后数值不正确,那么应该尝试用vodaconv_LE.exe(用的是Little Endian)重新生成raw文件。而一些可视化程序自身同时支持多种字节记录顺序,可在载入文件时选择。


5 vodaconv程序的原理

此程序源代码很容易读懂,没必要有太多注释。名为defvar的module用于保存格点设定。orgx,orgy,orgz是原点三个分量,dx,dy,dz是三个方向格点平移距离,nx,ny,nz是三个方向的点数。数组a保存原子信息,包括坐标,电荷,原子序数。cubmati1、cubmati2和cubmatr4分别是8bit含符号整型、16bit含符号整型和单精度浮点型动态分配空间的数组,用于储存格点数据。

下面介绍一下程序流程

1 载入体数据文件:首先根据输入文件的扩展名判定输入文件的格式。由于raw文件不含格点设定信息,而且数据类型未知,故需要由用户输入。平移矢量和原点坐标可输入可不输入,按ENTER可以使用默认值,这会将坐标原点设定到体数据的中心位置。.dx、.grd和.cub文件格式大同小异,载入.cub文件的代码在《Gaussian型cube文件简介及读、写方法和简单应用》(http://sobereva.com/125)一文中已经仔细介绍过了。这三种格式都是以浮点数表示数据,由于数据精度都不很高,为节约内存,用单精度浮点数记录就足够了。载入的数据暂存在cubmatr4数组中。

对于raw格式,若是单精度浮点型就直接载入到cubmatr4数组中。如果是8bit和16bit无符号整型,就分别载入到cubmati1和cubmati2中,然后将其数值转移到cubmatr4中,然后根据第三节所述对cubmatr4的负值数据部分进行运算以还原成原本的数值。另外,这也就是说在vodaconv程序内部对任何数据类型一律以单精度浮点方式保存,因为单精度浮点的储存精度和数值范围都远远高于8/16bit整型。

载入完毕后会显示体数据的格点设定信息以及最大和最小值。

2 对载入的数据的数值范围进行变换:这是给用户提供一个便利,可以避免数据越界,也可以将某一数值区域进行“放大”以使感兴趣的数值范围容易展现。在第一节也已经讨论了。如果想对数据做更丰富、灵活的变换操作,建议用寡人编写的Multiwfn程序2.3及以后版本(http://sobereva.com/multiwfn),它的主功能13专门用来操纵和提取cube格式的格点数据。

3 输出体数据文件:根据后缀名判断要输出的文件类型,将cubmatr4中的数据以相应方式输出即可。输出cube文件的代码在《Gaussian型cube文件简介及读、写方法和简单应用》里也已经详细讨论了。如本文第一节所述,可能会在cube文件中加入一个氢原子信息。