如何用NumPy处理Quantum Espresso能带数据且不丢失信息？

方法一：纯Python文本处理（无额外库依赖） #

直接按空行识别能带块边界，完整保留所有原始数据，无需预先知晓块行数：

import numpy as np
import matplotlib.pyplot as plt

def parse_qe_bands(filename):
    blocks = []
    current_block = []
    with open(filename, 'r') as f:
        for line in f:
            stripped_line = line.strip()
            # 跳过空行和注释行（QE能带文件开头可能有注释，可按需调整）
            if not stripped_line or stripped_line.startswith('#'):
                if current_block:
                    blocks.append(current_block)
                    current_block = []
                continue
            # 解析每行的k点和能量值
            k, energy = map(float, stripped_line.split())
            current_block.append([k, energy])
        # 处理文件末尾的最后一个能带块
        if current_block:
            blocks.append(current_block)
    
    # 提取第一块的完整k点坐标（保留重复值）
    k_points = np.array([row[0] for row in blocks[0]])
    # 提取每个能带的能量值，转置后适配绘图格式（每行对应一个k点，每列对应一个能带）
    band_energies = np.array([[row[1] for row in block] for block in blocks]).T
    
    return k_points, band_energies

# 使用示例
k, bands = parse_qe_bands('bands.dat')

# 绘制能带图
plt.figure(figsize=(8, 6))
for idx, band in enumerate(bands.T, 1):
    plt.plot(k, band, linewidth=1, label=f'Band {idx}')
plt.xlabel('k-path')
plt.ylabel('Energy (eV)')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()

核心优势 #

• 动态识别块边界：通过空行自动分割每个能带，完全适配不同系统的k点采样数量
• 无数据丢失：完整保留第一块的所有k点（包括重复的高对称点），不会像numpy.unique那样丢失信息
• 兼容QE输出：自动过滤注释行和空行，适配标准QE能带文件格式

方法二：用Pandas简化处理（适合熟悉Pandas的用户） #

利用Pandas的文本处理能力，更简洁地完成分块和数据提取：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 读取文件并按空行分割成能带块
with open('bands.dat', 'r') as f:
    raw_blocks = f.read().split('\n\n')

# 过滤无效块，将每个块转为DataFrame
band_dfs = []
for block in raw_blocks:
    lines = [line.strip() for line in block.split('\n') if line.strip() and not line.startswith('#')]
    if not lines:
        continue
    df = pd.DataFrame([map(float, line.split()) for line in lines], columns=['k', 'energy'])
    band_dfs.append(df)

# 提取k点和所有能带的能量矩阵
k_points = band_dfs[0]['k'].values
band_energies = np.column_stack([df['energy'].values for df in band_dfs])

# 绘图
plt.plot(k_points, band_energies, linewidth=1)
plt.xlabel('k-path')
plt.ylabel('Energy (eV)')
plt.show()

避坑提醒 #

• 禁止用numpy.unique处理k点：QE在高对称点（如Γ点）可能重复输出相同k坐标，对应不同能量值，去重会丢失关键数据
• 禁止固定行数提取：不同系统的k点采样数量差异很大，固定行数会导致代码完全无法通用，动态分块是唯一可靠方案

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