gromacs中xtc文件格式简介

本篇文章主要用于说明gromacs中xtc文件的编码格式。

文件整体构成

xtc文件一帧接着一帧，没有文件头、也没有末尾标记。读到文件结束即读完。每一帧的布局如下(单位：字节)：

字段	类型	大小	说明
magic	int32	4	魔数，有1995和2023两种可取的值
natoms	int32	4	本帧原子的个数，正常情况下全文件不变
step	int32	4	帧号(也就是步数)
time	float32	4	时间(ps)
box[9]	float32	36	盒子矩阵(3x3，按行存储)
坐标数据	变长	变长	随原子数决定格式

魔数字段

xtc头部的第一个字段就是魔数字段，该字段用来区分不同的版本。对应该字段的定义详见gromacs源代码中的src/gromacs/fileio/xdrf.h:52。

魔数	十进制(十六进制)	含义
XTC_MAGIC	1995(0x7cb)	传统格式，压缩缓冲区大小使用32位整数进行记录
XTC_NEW_MAGIC	2023(0x7e7)	新格式(用于原子数超过约3亿的系统)，压缩缓冲区大小使用64位整数记录

旧格式受限于32位整数，最多能存约2.98亿个原子(XTC_1995_MAX_NATOMS=298261617，详见gromacs源代码中的src/gromacs/fileio/xdrf.h:74)。当原子数超过该阈值时，gromacs自动使用新格式。

盒子矩阵字段

该字段按照行优先顺序依次为9个float，也就是存储顺序为box[0][0]、box[0][1]、box[0][2]、box[1][0]、box[1][1]、box[1][2]、box[2][0]、box[2][1]、box[2][2]。

坐标数据字段

坐标数据字段根据原子数的不同，可以分为两种情况。

情况1：原子数≤9

当natoms≤9时，坐标区的布局如下：

顺序	字段	编码	大小
1	lsize	int32	4字节
2	coordinate[3*natoms]	float32	12*natoms字节

也就是在此种情况下，坐标不进行压缩，直接写出坐标对应的IEEE 754单精度浮点数。lsize字段的取值为体系中原子的个数，也就是natoms的值。

情况2：原子数＞9

在此种情况下，坐标会进行压缩。坐标区的布局如下：

顺序	字段	编码	大小
1	lsize	int32	4字节
2	precision	float32	4字节
3	minint[0..2]	int32	12字节
4	maxint[0..2]	int32	12字节
5	smallidx	int32	4字节
6	buffer_size	int32(魔数为1995)或int64(魔数为2023)	4字节或8字节
7	buffer[]	原始字节	见下

上面的字段中，lsize与帧头里的natoms相同(写入时重复一次，用于校验)。buffer_size是压缩后位缓冲区的字节数。1995格式使用4字节，而2023格式使用8字节。buffer[]是通过xdr_opaque写入的。xdr_opaque会把字节数补足到4的整数倍，不足部分补0。也就是说，读入buffer_size个字节后，还要跳过(4-buffer_size%4)%4个填充字节才到达下一帧的头部。

gromacs中坐标压缩算法可以分为三个部分： - 定点化：把浮点坐标乘以精度四色五人成整数 - 范围编码：记录每维的最大值与最小值，把坐标转成相对于最小值的偏移，以此压缩位数 - 差分+游程编码：相邻原子坐标差很小时，只存差值，并动态调整差值的位数

定点化

对于每一个坐标分量x:

lf = x * precision;
if (lf >= 0) lf += 0.5; else lf -= 0.5;  // 四舍五入(正数加0.5，负数减0.5)
i = (int)lf;

范围编码

对三个维度分别统计整数坐标的最小值minint[d]与最大值maxint[d]，然后：

sizeint[d] = maxint[d] - minint[d] + 1; // 该维需要的取值个数

把坐标转换为相对偏移：tmpcoord[d] = coord_int[d] - minint[d]，取值范围是[0,sizeint[d])。需要注意的是，minint[d]是有符号int32。

根据sizeint[d]值的不同，编码相对偏移有两种模式： - 模式1：小范围。如果三个sizeint[d]都<=0xFFFFFF(2^24-1)。那么把它们打包成一个多字节的整数。三个偏移在同一个大整数V中按x高位、y中位以及z低位排布。

若a0、a1和a2分别为x、y和z三个方向的相对偏移，则

V = a0 * sizeint[1] * sizeint[2] + a1 * sizeint[2] + a2

读取时，可使用下面的方法恢复：

a2 = V % sizeint[2]
a1 = (V / sizeint[2]) % sizeint[1]
a0 = V / (sizeint[1] * sizeint[2])

- 模式2：大范围。如果任一sizeint[d]>0xFFFFFF，则每一维单独写入，不打包成多字节整数。

差分+游程编码

坐标编码的核心技巧：分子里相邻原子坐标往往很接近，存储当前原子相对于前一个原子的差值比存储绝对值省很多位。gromacs中利用magicints[]这张表来动态管理小差值的位数。

magicints = {0,0,0,0,0,0,0,0,0, // 下标 0..8
8, 10, 12, 16, 20, 25, 32, 40, 50, 64, 80, // 下标 9..19
101, 128, 161, 203, 256, 322, 406, 512, 645, 812,
1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501, 8192,
10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536, 82570,
104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287,
660561, 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304,
5284491, 6658042, 8388607, 10568983, 13316085, 16777216} // 下标 72

上面的magicints表的表长为73，第一个非0下标为9。下标0~8全为0,不使用。magicints[]本质上是把位数预算和分量取值范围互相换算的查找表。任何三个<magicints[i]的分量打包后都能用i位装下。

下面是关于差分+游程编码的相关字段。

字段	位置/位数	含义与作用
flag	每个大坐标之后，1位	本原子的游程/档位是否发生变化的标记，若为1，发生变化，后面会紧跟着5位run5。若为0，不发生变化，沿用上一次的run
run5	仅在flag=1时存在，5位	把“游程长度”和“档位变化”打包进一个5位数
run	从run5恢复	游程长度，单位是分量(不是原子数)
is_smaller	从run5恢复	档位调整量，取值为-1或0或1。1:升档(小差值用更多位、阈值smallnum变大)、-1:降档、0:不变
smallidx	帧头写入(int32)	magicints表的下标，数值上恰好等于“每个小差值分量所占的位数”

案例解析

情况1：原子数≤9

在这里使用spc2-traj.xtc这个含有6个原子的xtc文件进行解析。下面是该文件中内容的十六进制表示：

----------- 第一帧 --------------
00 00 07 cb      → magic = 0x7cb = 1995
00 00 00 06      → natoms = 6
00 00 00 00      → step = 0
00 00 00 00      → time = 0.0
40 40 a3 d7      → box[0][0] = 0x4040a3d7 = 3.00999
00 00 00 00      → box[0][1] = 0
00 00 00 00      → box[0][2] = 0
00 00 00 00      → box[1][0] = 0!
40 40 a3 d7      → box[1][1] = 0x4040a3d7 = 3.00999
00 00 00 00      → box[1][2] = 0
00 00 00 00      → box[2][0] = 0
00 00 00 00      → box[2][1] = 0
40 40 a3 d7      → box[2][2] = 0x4040a3d7 = 3.00999
00 00 00 06      → lsize = 6
3f 11 a9 fc      → coordinate1_x = 0x3f11a9fc = 0.56900
3f a3 33 33      → coordinate1_y = 0x3fa33333 = 1.27499
3f 95 1e b8      → coordinate1_z = 0x3f951eb8 = 1.16499
3e f3 b6 46      → coordinate2_x = 0x3ef3b646 = 0.47600
3f a2 4d d3      → coordinate2_y = 0x3fa24dd3 = 1.26800
3f 90 62 4e      → coordinate2_z = 0x3f90624e = 1.12800
...
...
3f bf 7c ee      → coordinate6_x = 0x3fbf7cee = 1.49600
3f c2 b0 21      → coordinate6_y = 0x3fc2b021 = 1.52100
3f 1f 7c ee      → coordinate6_z = 0x3f1f7cee = 0.62300
----------- 第二帧 --------------
00 00 07 cb      → magic = 0x7cb = 1995
00 00 00 06      → natoms = 6
00 00 00 01      → step = 1 
3f 80 00 00      → time = 1.0
...
...

情况2：原子数＞9

在这里使用spc2-traj.xtc这个含有6个原子的xtc文件进行解析。下面是该文件中内容的十六进制表示(//后的为注释)：

----------- 第一帧 --------------
00 00 07 cb      → magic = 0x7cb = 1995
00 00 00 9c      → natoms = 156
00 00 00 00      → step = 0
00 00 00 00      → time = 0.0
40 bc ff 97      → box[0][0] = 0x40bcff97 = 5.90619
00 00 00 00      → box[0][1] = 0
00 00 00 00      → box[0][2] = 0
00 00 00 00      → box[1][0] = 0
40 db 0b 0f      → box[1][1] = 0x40db0b0f = 6.84509
00 00 00 00      → box[1][2] = 0
00 00 00 00      → box[2][0] = 0
00 00 00 00      → box[2][1] = 0
40 43 4f 0e      → box[2][2] = 0x40434f0e = 3.05170
00 00 00 9c      → lsize = 156
44 7a 00 00      → precision = 1000.0
00 00 0b f4      → minint_x = 3060
00 00 02 72      → minint_y = 626
00 00 00 0a      → minint_z = 10
00 00 10 40      → maxint_x = 4160
00 00 08 f0      → maxint_y = 2288
00 00 0b e2      → maxint_z = 3042
// x维需要的取值个数：4160 - 3060 + 1 = 1101
// y维需要的取值个数：2288 - 626 + 1 = 1663
// z维需要的取值个数：3042 - 10 + 1 = 3033
// 那么，有1101*1663*3033=5,553,310,779种情况，那么需要33位进行存储

00 00 00 15      → smallidx = 21
00 00 02 48      → buffer_size = 584 (584 % 4 == 0，无填充)

// 下面为buffer[]
c0 10 65 8f  
// 组装成一个大整数V = 0xc0|0x10<<8|0x65<<16|0x8f<<24|0<<32 
// 也就是1000 1111 0110 0101 0001 0000 1100 0000 = 2,405,765,312
// a2 = 2405765312 % 3033 = 1844
// a1 = (2405765312 / 3033) % 1663 = 793196 % 1663 = 1608
// a0 = 2405765312 / (3033 * 1663 ) = 476
// 加上minint得到整数坐标(3536,2234,1854)
// 除以precision=1000得(3.536, 2.234, 1.854)

43 e9 92 81  // 随后，读取一位flag位，为1。若flag为1，再读5位run5
4c 89 da 32  // atom0的原子坐标对应的位为0xf4、0xc9、0x40、0xa6，最后一位为0
// 0xf4 | 0xc9<<8 | 0x40<<16 | 0xa6<<24 = 0xa640c9f4 = 2,789,263,860
// a2 = 2789263860 % 3033 = 1806
// a1 = (2789263860 / 3033) % 1663 = 1662
// a0 = 2789263860 / (3033 * 1663) = 552
// 加上minint得到整数坐标(3612, 2288, 1816)
// 除以precision=1000得(3.612, 2.288, 1.816)

整理可以得到如下的表：

原子	大坐标位区间	偏移量	整数坐标	flag的位次	flag的值	run5的位次	run5的值	run	is_smaller	smallidx
0	[0, 33)	(476, 1608, 1844)	(3536, 2234, 1854)	33	1	[34, 38]	1	0	0	21→21
1	[39, 72)	(552, 1662, 1806)	(3612, 2288, 1816)	72	1	[73, 77]	2	0	+1	21→22
2	[78, 111)	(410, 1558, 1772)	(3470, 2214, 1782)	111	1	[112, 116]	2	0	+1	22→23

其中，在写端，run5 = run + is_smaller + 1，在读端，is_smaller = run5 % 3; run -= is_smaller; is_smaller -= 1。

评论