读取Gromacs的.gro大文件写入.h5格式报错的解决方案求助

错误根因说明 #

• ValueError（列数不匹配）：.gro是Gromacs的分帧轨迹格式，每帧有固定结构，不是所有行都是6列的原子数据：
  1. 每帧第1行：帧标题（含时间信息，拆分后列数不定，你报错的7列行就是该类行）
  2. 每帧第2行：当前帧原子总数（拆分后仅1列，对应报错的1列行）
  3. 接下来N行（N为原子总数）：原子数据行，固定6列
  4. 每帧最后1行：模拟盒子矢量参数，拆分后为3/9列
     你之前无差别拆分所有行按原子数据处理，自然会触发列数不匹配错误。
• TypeError（字符串无法写入H5）：一是拆分后得到的字符串未按字段类型转换（原子数据前2列为字符串、第3列为整数、后3列为浮点数）；二是h5py默认不支持直接写入numpy变长字符串类型（<U*），需要显式指定固定长度的字符串dtype。
• 额外致命问题：用readlines()读取50G文件会直接占满系统内存，完全无法运行；且循环结束后你仅保留了最后一行的拆分结果，所有历史数据全部丢失，写入H5的内容完全不符合预期。

最优读取方案（推荐） #

直接使用专门处理分子动力学轨迹的MDAnalysis库，内置gro格式解析能力，支持流式读取大轨迹文件，不会触发内存溢出，也不需要自己处理格式解析的坑。

实现代码 #

import h5py
import numpy as np
import MDAnalysis as mda

# 加载gro轨迹文件，流式读取不会加载全量数据到内存
u = mda.Universe("pep.gro", in_memory=False)
n_frames = len(u.trajectory)
n_atoms = u.atoms.n_atoms

# 初始化H5文件，创建可扩展数据集
with h5py.File("pep1.h5", "w") as hdf:
    # 存储原子元信息：残基名、原子名、原子编号，仅存一次即可
    hdf.create_dataset("resname", data=np.array(u.atoms.resnames, dtype="S10"))
    hdf.create_dataset("atomname", data=np.array(u.atoms.names, dtype="S10"))
    hdf.create_dataset("atomid", data=u.atoms.ids.astype(np.int32))
    # 存储坐标数据集，按帧存储，shape为(帧数, 原子数, 3)
    coords_dset = hdf.create_dataset(
        "coords", 
        shape=(0, n_atoms, 3), 
        maxshape=(None, n_atoms, 3), 
        dtype=np.float32,
        chunks=True
    )
    # 存储每帧的时间
    time_dset = hdf.create_dataset(
        "time",
        shape=(0,),
        maxshape=(None,),
        dtype=np.float32
    )
    
    # 逐帧迭代读取，不会爆内存
    for i, ts in enumerate(u.trajectory):
        coords_dset.resize(i+1, axis=0)
        coords_dset[i] = ts.positions
        time_dset.resize(i+1, axis=0)
        time_dset[i] = ts.time

手动解析方案（无第三方依赖） #

如果不想安装MDAnalysis，可以自己按gro格式规则逐行流式解析：

import h5py
import numpy as np

# 提前确认原子数（可以先读前两行获取，避免后面重复判断）
with open("pep.gro", "r") as f:
    f.readline()
    n_atoms = int(f.readline().strip())

# 初始化H5数据集
with h5py.File("pep1.h5", "w") as hdf, open("pep.gro", "r") as f:
    coords_dset = hdf.create_dataset(
        "coords",
        shape=(0, n_atoms, 3),
        maxshape=(None, n_atoms, 3),
        dtype=np.float32,
        chunks=True
    )
    time_dset = hdf.create_dataset(
        "time",
        shape=(0,),
        maxshape=(None,),
        dtype=np.float32
    )
    frame_idx = 0
    coords_buffer = np.zeros((n_atoms, 3), dtype=np.float32)
    
    while True:
        # 读标题行，提取时间
        title_line = f.readline()
        if not title_line:
            break # 文件读取结束
        time = float(title_line.strip().split("t=")[-1])
        # 读原子数行，跳过
        f.readline()
        # 读n_atoms行原子数据
        for i in range(n_atoms):
            parts = f.readline().strip().split()
            coords_buffer[i] = list(map(float, parts[3:6]))
        # 读盒子行，跳过
        f.readline()
        # 写入H5
        coords_dset.resize(frame_idx+1, axis=0)
        coords_dset[frame_idx] = coords_buffer
        time_dset.resize(frame_idx+1, axis=0)
        time_dset[frame_idx] = time
        frame_idx += 1

内容的提问来源于stack exchange，提问作者Mahesh