使用 AbiPy 提取 ABINIT 几何优化轨迹并在 OVITO 中可视化

痛点及场景动机分析

在材料结构优化与分子动力学模拟中，ABINIT 生成的结果与过程主要记录在 HIST.nc 文件内。尽管研究者可以直接利用标准库读取静态的原子坐标，但在复原多步优化、含有晶胞周期性波动的完整演化过程并直观呈现时，传统的解析手段往往无法做到直接衔接。

此时，OVITO 作为一款出色的三维原子级可视化软件显现出了绝佳适配性。然而，为了避免因边界信息缺失导致跨晶格结构的视觉错误，喂给 OVITO 的数据必须具有准确的晶胞（Lattice）规范。在此背景下，利用 pymatgen 以及 AbiPy 重新构造含特殊 Properties 头的 Extended XYZ 格式就成了一套必须而有效的工作流。

工具链工作流与核心代码实现解析

整个 Python 处理脚本涵盖了从运行依赖的热修复，到核心原子对象处理，最终统一序列化落盘等关键流程：

高版本 SciPy 的补丁兼容

在全新发布的 SciPy 版本中，计算积分的函数 cumtrapz 与 simps 已更名为 cumulative_trapezoid 与 simpson。为了防止 AbiPy 底层引用旧名而遭遇异常崩溃，我们在代码顶部嵌入了动态重定向功能予以快速修复兼容：

# 补丁：适配高版本 scipy
import scipy.integrate as _sci_int
if not hasattr(_sci_int, 'cumtrapz'):
    _sci_int.cumtrapz = _sci_int.cumulative_trapezoid
if not hasattr(_sci_int, 'simps'):
    _sci_int.simps = _sci_int.simpson

构建并序列化 Extended XYZ 数据元

为了适应 OVITO 原生协议的加载要求，系统将基于 pymatgen.Structure 参数化抽取对应的 $3 \times 3$ 晶格向量，使用行优先平铺法则整理为一段向量标量，并附着到代表当前步的 Lattice 及 Properties 首部信息条目之上：

def write_extxyz(structures, filepath):
    """
    将 pymatgen.Structure 列表转换为 Extended XYZ 格式。
    Extended XYZ 格式包含 Lattice 关键字，OVITO 能够自动识别并绘制正确的晶胞/真空层盒子。
    """
    with open(filepath, 'w') as f:
        for step, struct in enumerate(structures):
            f.write(f"{len(struct)}\n")

            # 提取 3x3 晶格向量，按行展平
            lat = struct.lattice.matrix
            lat_str = f"{lat[0,0]} {lat[0,1]} {lat[0,2]} {lat[1,0]} {lat[1,1]} {lat[1,2]} {lat[2,0]} {lat[2,1]} {lat[2,2]}"

            # 使用 OVITO 支持的 Properties 和 Lattice 关键字
            info_line = f'Lattice="{lat_str}" Properties=species:S:1:pos:R:3 Step={step}'
            f.write(info_line + '\n')

            # 写入原子和对应的笛卡尔坐标
            for site in struct:
                f.write(f"{site.specie.symbol} {site.x:.6f} {site.y:.6f} {site.z:.6f}\n")

批量提取与自动化控制

利用命令行界面接收指定工作路径，结合模式化搜索以捕获多文件目录树下面所有的 .nc 源计算归档，随后程序可以一键驱动并实时反馈进度情况与错误追踪：

# 批量读取并导出逻辑片断
hist_files = sorted(glob.glob(os.path.join(args.calc_dir, '*', '*HIST.nc')))
for hist_path in hist_files:
    dir_name = os.path.basename(os.path.dirname(hist_path))
    out_path = os.path.join(args.out_dir, f"{dir_name}_traj.extxyz")

    try:
        hist = abilab.abiopen(hist_path)
        structures = hist.structures
        write_extxyz(structures, out_path)
        hist.close()
    except Exception as e:
        print(f"解析错误: {e}")

完整代码合集分享

以下为完整可复用脚本 extract-opt-traj-abipy.py 的代码分享，请在运行中保障包含了 AbiPy、Pymatgen 以及匹配计算工具的环境状态（如加载虚拟 conda dft 环境）：

#!/usr/bin/env python3
import os
import glob
import argparse

# 补丁：适配高版本 scipy
import scipy.integrate as _sci_int
if not hasattr(_sci_int, 'cumtrapz'):
    _sci_int.cumtrapz = _sci_int.cumulative_trapezoid
if not hasattr(_sci_int, 'simps'):
    _sci_int.simps = _sci_int.simpson

# 导入 abipy
try:
    import abipy.abilab as abilab
except ImportError:
    raise ImportError("Failed to import abipy. Activate conda environment 'dft' first.")

def write_extxyz(structures, filepath):
    """
    将 pymatgen.Structure 列表转换为 Extended XYZ 格式。
    Extended XYZ 格式包含 Lattice 关键字，OVITO 能够自动识别并绘制正确的晶胞/真空层盒子。
    """
    with open(filepath, 'w') as f:
        for step, struct in enumerate(structures):
            f.write(f"{len(struct)}\n")

            # 提取 3x3 晶格向量，按行展平
            lat = struct.lattice.matrix
            lat_str = f"{lat[0,0]} {lat[0,1]} {lat[0,2]} {lat[1,0]} {lat[1,1]} {lat[1,2]} {lat[2,0]} {lat[2,1]} {lat[2,2]}"

            # 使用 OVITO 支持的 Properties 和 Lattice 关键字
            info_line = f'Lattice="{lat_str}" Properties=species:S:1:pos:R:3 Step={step}'
            f.write(info_line + '\n')

            # 写入原子和对应的笛卡尔坐标
            for site in struct:
                f.write(f"{site.specie.symbol} {site.x:.6f} {site.y:.6f} {site.z:.6f}\n")

def main():
    parser = argparse.ArgumentParser(description="使用 abipy 提取几何优化结构，包含晶胞信息")
    parser.add_argument('--calc-dir', default='calc', help='计算目录')
    parser.add_argument('--out-dir', default='str/opt_traj', help='轨迹输出目录')
    args = parser.parse_args()

    os.makedirs(args.out_dir, exist_ok=True)

    hist_files = sorted(glob.glob(os.path.join(args.calc_dir, '*', '*HIST.nc')))

    if not hist_files:
        print(f"No HIST.nc found in {args.calc_dir}")
        return

    n_ok = 0
    n_fail = 0

    for hist_path in hist_files:
        dir_name = os.path.basename(os.path.dirname(hist_path))
        # 使用 .extxyz 后缀（或 .xyz），并在文件头加上盒子参数
        out_path = os.path.join(args.out_dir, f"{dir_name}_traj.extxyz")

        try:
            hist = abilab.abiopen(hist_path)
            structures = hist.structures
            write_extxyz(structures, out_path)
            hist.close()
            print(f"[OK  ] {dir_name:<25} 提取了 {len(structures):>4} 步 -> {out_path}")
            n_ok += 1
        except Exception as e:
            print(f"[FAIL] {dir_name:<25} Error: {e}")
            n_fail += 1

    print(f"\nDone! Success: {n_ok}, Fail: {n_fail}. Saved to {args.out_dir}/")

if __name__ == '__main__':
    main()