求助：如何用Python绘制含PME数据的DX文件平均势能剖面

依赖安装 #

先安装必要的Python库：

pip install numpy matplotlib

完整代码 #

import numpy as np
import matplotlib.pyplot as plt

def parse_dx_file(file_path):
    """解析DX格式的势能文件，返回网格参数和三维势能数组"""
    with open(file_path, 'r') as f:
        lines = [line.strip() for line in f if line.strip() and not line.startswith('#')]
    
    # 解析网格维度
    counts_line = [line for line in lines if 'gridpositions counts' in line][0]
    nx, ny, nz = map(int, counts_line.split()[-3:])
    
    # 解析原点坐标
    origin_line = [line for line in lines if 'origin' in line][0]
    x0, y0, z0 = map(float, origin_line.split()[-3:])
    
    # 解析网格步长
    delta_lines = [line for line in lines if 'delta' in line]
    dx = float(delta_lines[0].split()[-3])
    dy = float(delta_lines[1].split()[-2])
    dz = float(delta_lines[2].split()[-1])
    
    # 找到数据起始位置
    data_start_idx = next(i for i, line in enumerate(lines) if 'data follows' in line) + 1
    data_lines = lines[data_start_idx:]
    
    # 读取并整理势能数据
    potential_data = []
    for line in data_lines:
        potential_data.extend(map(float, line.split()))
    potential_3d = np.array(potential_data).reshape(nz, ny, nx)  # 若结果异常，可尝试调整维度顺序
    
    # 生成坐标轴
    x = np.linspace(x0, x0 + (nx-1)*dx, nx)
    y = np.linspace(y0, y0 + (ny-1)*dy, ny)
    z = np.linspace(z0, z0 + (nz-1)*dz, nz)
    
    return x, y, z, potential_3d

def calculate_average_profile(potential_3d, axis='z'):
    """沿着指定轴计算平均势能剖面"""
    if axis == 'x':
        avg_profile = np.mean(potential_3d, axis=(0,1))  # 对y、z轴取平均
    elif axis == 'y':
        avg_profile = np.mean(potential_3d, axis=(0,2))  # 对x、z轴取平均
    elif axis == 'z':
        avg_profile = np.mean(potential_3d, axis=(1,2))  # 对x、y轴取平均
    else:
        raise ValueError("axis must be 'x', 'y', or 'z'")
    return avg_profile

# 主流程
if __name__ == '__main__':
    # 替换为你的DX文件路径
    dx_file_path = 'your_pme_potential.dx'
    
    # 解析DX文件
    x, y, z, potential_3d = parse_dx_file(dx_file_path)
    
    # 计算指定轴的平均势能剖面（可替换为'x'或'y'）
    avg_profile = calculate_average_profile(potential_3d, axis='z')
    
    # 绘制剖面
    plt.figure(figsize=(8, 5))
    plt.plot(z, avg_profile, linewidth=2, color='b')
    plt.xlabel('Z Coordinate (Å)')  # 根据数据实际单位调整，比如nm
    plt.ylabel('Average Potential (kJ/mol)')  # 单位按需修改，比如kcal/mol
    plt.title('PME Average Potential Profile along Z-axis')
    plt.grid(alpha=0.3)
    plt.tight_layout()
    plt.show()

注意事项 #

• 维度顺序：DX文件的存储顺序可能存在差异，如果绘图结果异常，可尝试将reshape(nz, ny, nx)改为reshape(nx, ny, nz)或其他组合。
• 单位适配：代码中的坐标和势能单位是常见模拟场景的默认值，请根据你的实际模拟数据修改标签内容。
• 格式兼容：如果你的DX文件包含特殊注释或非标准格式，需要微调parse_dx_file函数中的筛选逻辑，比如修改行匹配的关键词。

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