从gromacs文件中读取数据并写入hdf5文件格式

Question

我正在尝试逐行从 .gro 文件中读取数据，并希望将其写入 .h5 文件格式的数据。但是出现类型错误："No conversion path ford type: type('<U7')"。我猜读取的数据是字符串格式。我尝试使用 np.arrays 将其转换为数组，但它不起作用。谁能帮我解决这个问题？或者有更好的方法来读取数据吗？我无法使用 np.loadtxt，因为数据大小约为 50GB。

.gro 文件的格式如下所示

Generated by trjconv : P/L=1/400 t=   0.00000
11214
    1P1     aP1    1  80.48  35.36   4.25
    2P1     aP1    2  37.45   3.92   3.96
Generated by trjconv : P/L=1/400 t=   10.00000
11214
    1P1     aP1    1  80.48  35.36   4.25
    2P1     aP1    2  37.45   3.92   3.96
Generated by trjconv : P/L=1/400 t=   20.00000
11214
    1P1     aP1    1  80.48  35.36   4.25
    2P1     aP1    2  37.45   3.92   3.96
Generated by trjconv : P/L=1/400 t=   30.00000
11214
    1P1     aP1    1  80.48  35.36   4.25
    2P1     aP1    2  37.45   3.92   3.96
Generated by trjconv : P/L=1/400 t=   40.00000
11214
    1P1     aP1    1  80.48  35.36   4.25
    2P1     aP1    2  37.45   3.92   3.96

错误：

ValueError: Some errors were detected !
    Line #5 (got 7 columns instead of 6)
    Line #6 (got 1 columns instead of 6)
    Line #9 (got 7 columns instead of 6)
    Line #10 (got 1 columns instead of 6)
    Line #13 (got 7 columns instead of 6)
    Line #14 (got 1 columns instead of 6)
    Line #17 (got 7 columns instead of 6)
    Line #18 (got 1 columns instead of 6)

这是我的小代码：

import h5py
import numpy as np
# First step is to read .gro file
f = open('pep.gro', 'r')
data = f.readlines()
for line in data:
    reading = line.split()
    #print(type(reading))
    #dat = np.array(reading).astype(int)

# Next step is to write the data to .h5 file
with h5py.File('pep1.h5', 'w') as hdf:
    hdf.create_dataset('dataset1', data=reading)

Answer 1

您首先创建具有大量行 [shape=(1_000_000)] 的 HDF5 数据集，然后使用 maxshape 参数使其可扩展。 maxshape=(None,) 的值将允许无限行。我定义了一个简单的 dtype 来匹配您的数据。如果需要，可以自动为不同的文件格式创建匹配的 dtype。

您遇到了 Unicode 错误，因为 h5py 不支持将字符串作为 Unicode 数据。（NumPy 默认从字符串创建 Unicode 数据。）解决此限制的方法是预先为数组定义一个数据类型（使用 'S#'，其中 NumPy 具有“

接下来使用 np.genfromtxt 将您的内容直接读入 NumPy 数组。使用 skip_header 和 max_rows 参数增量读取。包括 dtype 参数和用于创建上述数据集的数据类型。

为了测试增量读取，我将您的文件扩展到 54 行（用于 3 个读取循环）。出于性能原因，您需要使用更大的值来读取 50GB（将 incr 设置为您可以读入内存的内容——从 100_000 行开始）。

下面的代码：（修改为跳过前两行

import h5py
import numpy as np

#define a np.dtype for gro array/dataset (hard-coded for now)
gro_dt = np.dtype([('col1', 'S4'), ('col2', 'S4'), ('col3', int), 
                   ('col4', float), ('col5', float), ('col6', float)])

# Next, create an empty .h5 file with the dtype
with h5py.File('pep1.h5', 'w') as hdf:
    ds= hdf.create_dataset('dataset1', dtype=gro_dt, shape=(20,), maxshape=(None,)) 

    # Next read line 1 of .gro file
    f = open('pep.gro', 'r')
    data = f.readlines()
    ds.attrs["Source"]=data[0]
    f.close()

    # loop to read rows from 2 until end
    skip, incr, row0 = 2, 20, 0 
    read_gro = True
    while read_gro:
        arr = np.genfromtxt('pep.gro', skip_header=skip, max_rows=incr, dtype=gro_dt)
        rows = arr.shape[0]
        if rows == 0:
            read_gro = False 
        else:    
            if row0+rows > ds.shape[0] :
                ds.resize((row0+rows,))
            ds[row0:row0+rows] = arr
            skip += rows
            row0 += rows

Answer 2

读取2000个时间步长的过程类似，只是步长顺序不同。此外，由于我们知道每个时间步长的数据行数（从第二行 header 开始），这简化了数据集的创建。我们可以使用它来调整每个数据集的大小，并避免创建必须调整大小的可扩展数据集。

我正在写一个新的答案来说明如何做到这一点。关于 dtype 和 Numpy Unicode 字符串的评论仍然适用。此代码适用于我创建的具有 2 个时间步长的测试文件。修改 n_steps= 以使用 2000 步。

import h5py
import numpy as np

n_steps = 2 # 2000
csv_file = 'pep2.gro'

#define a np.dtype for gro array/dataset (hard-coded for now)
gro_dt = np.dtype([('col1', 'S4'), ('col2', 'S4'), ('col3', int), 
                   ('col4', float), ('col5', float), ('col6', float)])

# Open the file for reading and
# create an empty .h5 file with the dtype above   
with open(csv_file, 'r') as f, \
     h5py.File('pep2.h5', 'w') as hdf:

    data = f.readlines()
    skip= 2

    for step in range(n_steps):
        # Read text header line for THIS time step
        header = data[skip-2]
        # get number of data rows
        no_rows = int(data[skip-1]) 
        arr = np.genfromtxt(csv_file, skip_header=skip, max_rows=no_rows, dtype=gro_dt)
        if arr.shape[0] > 0:
            # create a dataset for THIS time step
            ds= hdf.create_dataset(f'dataset_{step:04}', data=arr) 
            #create attributes for this dataset / time step
            hdr_tokens = header.split()
            ds.attrs['raw_header'] = header
            ds.attrs['Generated by'] = hdr_tokens[2]
            ds.attrs['P/L'] = hdr_tokens[4].split('=')[1]
            ds.attrs['Time'] = hdr_tokens[6]
        skip += no_rows + 3

从gromacs文件中读取数据并写入hdf5文件格式

Reading data from gromacs file and write it to the hdf5 file format

python

hdf5

h5py