用于自动提取芯轨道(分子和原子轨道)的程序

wangsy 发布于 2025-12-5 14:47 阅读：47 密度泛函

`extract_orbital_composition.py` 程序使用说明

本程序由wangsy于2025年11月19日完成。
如有需要，请联系wang.sy@vip.163.com 。

1. 简介 (Program Overview)

extract_orbital_composition.py 是一个用于量子化学计算的后处理工具。它旨在自动化地从.fchk (formatted checkpoint) 文件中提取指定原子在占据分子轨道中的成分百分比。

该脚本的核心功能是调用业界知名的波函数分析软件 Multiwfn，向其发送一系列自动化指令，然后解析其文本输出，最后将结果以清晰的表格或缓存文件的形式呈现给用户。

关键特性

自动化 Multiwfn 调用: 自动运行 Multiwfn 并执行轨道成分分析（菜单 8 -> 8 -> 1 -> -2）。
计算类型自动检测: 自动从 .fchk 文件中读取电子数和轨道信息，以判断计算是限制性 (Restricted, RHF/ROHF) 还是非限制性 (Unrestricted, UHF)。
MOKIT 自动局域化: 当脚本检测到轨道高度离域时（即，目标原子的最高轨道贡献低于 95%），它能自动调用 MOKIT 工具包执行 Pipek-Mezey (PM) 轨道局域化，并使用局域化后的新 fchk 文件重新进行分析。
灵活的输出:
1. 默认在屏幕上打印格式化的表格。
2. 可使用 -o 将表格保存到文件。
3. 可使用 -c 生成用于 qchem_xray.py 脚本的专用缓存文件。

⚠️ 重要：运行依赖

在运行此脚本之前，您必须确保满足以下环境要求：

Multiwfn: Multiwfn 的可执行文件必须位于您系统的 PATH 环境变量中。脚本会通过在 shell 中执行 Multiwfn 命令来调用它。
MOKIT (用于局域化):
- 您必须有一个名为 mokit-py39 的 conda 环境。
- 此环境中必须已正确安装 mokit 库及其依赖项。
- 脚本通过 conda run -n mokit-py39 ... 命令来执行局域化相关的 Python 和 fch2qchem 指令。

2. 逻辑和流程 (Logic and Flow)

该脚本的执行流程可以分为以下几个关键步骤：

参数解析: 脚本启动，读取用户提供的命令行参数（例如 fchk_file, atom_number, threshold 等）。
环境检查: 脚本首先检查 Multiwfn 是否在 PATH 中。如果找不到，程序将报错并退出。
文件初次分析 (Original FCHK):
- 脚本读取 .fchk 文件，提取 Alpha 电子数 (n_alpha)、Beta 电子数 (n_beta) 和总基函数数 (total_orbitals)。
- 通过检查文件中是否存在 "Beta MO coefficients" 等关键词，判断是限制性还是非限制性计算。
- 调用 Multiwfn: 脚本通过 subprocess 启动 Multiwfn [fchk_file] -isilent 进程。
- 发送指令: 它向 Multiwfn 的标准输入 (stdin) 发送一系列命令，以分析所有占据轨道 (RHF: 1 到 n_alpha；UHF: 1 到 n_alpha 以及 Beta 轨道 1 到 n_beta)。
- 解析输出: 脚本捕获 Multiwfn 的标准输出 (stdout)，并使用正则表达式解析表格，提取所有贡献超过指定阈值 (-t，默认 40%) 的轨道信息（编号、能量、占据数、百分比等）。
离域化检查:
- 脚本检查上一步中找到的所有轨道的最大贡献值 (max_comp)。
- 如果 0.0 < max_comp < 95.0% (95% 为 LOCALIZATION_TRIGGER_THRESHOLD)，则认为轨道高度离域，触发局域化流程。
MOKIT 局域化工作流 (可选):
- 如果触发了局域化（且用户未指定 --no-localize）：
- 交互式询问: 脚本会暂停并询问用户是否执行 MOKIT 局域化（y/n）。（如果使用 --force-localize 则跳过询问）。
- 执行局域化: 如果用户同意，脚本将：
  i. 生成一个临时的 Python 脚本，该脚本导入 mokit 库。
  ii. 调用 conda run -n mokit-py39 python ... 来执行 standardize_fch (标准化 fchk) 和 loc (PM 局域化)。
  iii. 这会生成一个新的 *_std_LMO.fch 文件。
  iv. 调用 conda run -n mokit-py39 fch2qchem ... 将新的 LMO fchk 转换为 Q-Chem save 目录。
  v. 将生成的 *_std_LMO 目录重命名为 save（如果 save 已存在，则备份为 save.bak）。
二次分析 (LMO FCHK):
- 如果成功执行了局域化，脚本将使用新生成的 *_std_LMO.fch 文件重复步骤 3。
- 此时，分析的是局域化分子轨道 (LMO)。
结果输出:
- 脚本将最终的分析结果（来自 LMO fchk 或原始 fchk）进行格式化。
- 标准模式: 打印一个详细的表格到屏幕（或用 -o 保存到文件）。
- 缓存模式 (-c): 生成一个 .cache 文件。
  - 标准: 2 列 (Alpha, Beta)。
  - 局域化后: 4 列 (LMO_Alpha, LMO_Beta, Orig_Alpha, Orig_Beta)，并包含特殊的页脚标记。

3. 程序的运行完整示例

假设我们有一个名为 complex.fchk 的非限制性 (UHF) 计算文件，我们关心的是 1 号原子（例如一个金属中心）的轨道贡献。

示例 1: 基本用法（打印到屏幕）

命令 (使用 10% 作为阈值):

python extract_orbital_composition.py complex.fchk 1 -t 10

输出 (无局域化):

如果轨道局域性良好（最高贡献 > 95%），脚本将直接打印分析结果。

==========================================================================================
                                原子 1 的轨道组成分析
==========================================================================================
文件: complex.fchk
计算类型: 非限制性 (Unrestricted)
Composition阈值: 10.0%
(注: 仅显示占据轨道)
==========================================================================================
注: 非限制性计算有两套独立的Alpha和Beta空间轨道。
------------------------------------------------------------------------------------------

Alpha 轨道:
------------------------------------------------------------------------------------------
      序号 类型                能量(a.u.)    占据     组成(%)         布居
------------------------------------------------------------------------------------------
       105 Alpha              -0.2512   1.000      98.501   1.000000
       100 Alpha              -0.4500   1.000      85.123   1.000000
        95 Alpha              -0.7800   1.000      15.300   1.000000

Beta 轨道:
------------------------------------------------------------------------------------------
      序号 类型                能量(a.u.)    占据     组成(%)         布居
------------------------------------------------------------------------------------------
    104(204) Beta               -0.2300   1.000      97.110   1.000000
     99(199) Beta               -0.4300   1.000      82.000   1.000000
------------------------------------------------------------------------------------------
注: Beta轨道序号格式为 真实序号(Multiwfn内部编号)。
------------------------------------------------------------------------------------------

示例 2: 交互式局域化

命令 (使用默认 40% 阈值):

python extract_orbital_composition.py complex.fchk 1

交互式提示:

假设脚本在第一次分析中发现最高贡献只有 85%（低于 95%）。

处理非限制性体系 (Unrestricted)...
 → 扫描范围 (Alpha,Beta): 1-105,201-204

警告: 排序第一的轨道百分比 (85.12%) 低于 95.0%。
      这可能表明轨道高度离域。

是否运行 MOKIT (PM) 局域化工作流? (y/n):

用户输入 y:

脚本将开始执行 MOKIT 工作流，并显示 MOKIT 的输出：

y
信息: 用户选择执行局域化。

============================================================
开始 MOKIT 局域化工作流 (目标原子: 1)
============================================================
--- [1/3] 正在 MOKIT 环境中运行 Python 脚本 (_temp_mokit_localize_...py)...
--- 正在导入 MOKIT...
--- [MOKIT 1/2] 正在标准化 complex.fchk -> complex_std.fch...
--- [MOKIT 2/2] 正在局域化 complex_std.fch (PM, 轨道 0-104)...
--- MOKIT Python 脚本执行完毕.
--- [2/3] 正在 MOKIT 环境中运行 fch2qchem (转换LMO fch为Q-Chem 'save' 目录)...
(fch2qchem output...)
--- [3/3] 正在重命名文件...
============================================================
局域化完成。生成新的fchk文件: complex_std_LMO.fch
============================================================

--- 正在分析局域化后的文件: complex_std_LMO.fch
处理非限制性体系 (Unrestricted)...
 → 扫描范围 (Alpha,Beta): 1-105,201-204

最终输出:

脚本现在打印的是对 complex_std_LMO.fch 文件的分析结果，轨道贡献现在应该更接近 100%。

==========================================================================================
                                原子 1 的轨道组成分析
==========================================================================================
文件: complex_std_LMO.fch
计算类型: 非限制性 (Unrestricted)
...
Alpha 轨道:
...
       105 Alpha              -0.2512   1.000      99.850   1.000000
       100 Alpha              -0.4500   1.000      99.100   1.000000
...

示例 3: 生成缓存文件（自动局域化）

命令 (使用 -c 生成缓存，使用 --force-localize 避免交互式询问):

python extract_orbital_composition.py complex.fchk 1 -c --force-localize

输出:

...
(MOKIT 局域化流程输出)
...

轨道缓存文件已成功生成: /path/to/project/1GS_orbital.cache

1GS_orbital.cache 文件内容:

由于执行了局域化，缓存文件将包含 4 列，并带有特殊的页脚。

n_alpha=105       n_beta=104
105|99.85        104|99.11        105|85.12        104|82.00
100|99.10        99|82.00         100|82.00        99|15.30

# --------------------------------------------------
# 警告: 原始fchk文件轨道百分比 < 95.0% (Max: 85.12%)
# LOCALIZED: LMO_A | LMO_B | ORIG_A | ORIG_B
# --------------------------------------------------

4. 高级选项 (Advanced Options)

以下是脚本支持的所有命令行参数的详细分类说明。

核心参数

fchk_file
- 必需。输入的 .fchk 文件路径。
atom_number
- 必需。要分析的目标原子的编号（从 1 开始）。

阈值与输出

-t, --threshold THRESHOLD
- 设置轨道组成贡献的百分比阈值。
- 示例: -t 10 (只显示贡献 > 10% 的轨道)。
- 默认: 40.0
-o, --output OUTPUT
- 将标准输出的表格重定向到指定的文本文件。

缓存模式 (qchem_xray.py)

-c, --cache
- 激活缓存模式。不打印标准表格，而是生成一个 qchem_xray.py 兼容的缓存文件。
- 默认文件名: [atom_number]GS_orbital.cache (例如: 1GS_orbital.cache)。
--cache-name CACHE_NAME
- 与 -c 配合使用，指定自定义的缓存文件名。

MOKIT 局域化控制

--force-localize
- 当检测到轨道离域时 (max comp \< 95%)，自动执行 MOKIT 局域化，不进行交互式询问。
- 适用于在其他脚本中调用此程序。
--no-localize
- 禁止 MOKIT 局域化。即使检测到离域，也只使用原始 fchk 文件进行分析。
--keep-temp-files
- 保留 MOKIT 局域化过程中产生的中间文件（例如 *_std.fch, *_std_LMO.in）。
- 默认情况下，这些文件会被自动清理。此选项用于调试 MOKIT 流程。

自旋控制 (Spin Control)

-a, --alpha-only
- 仅分析和输出 Alpha 轨道。
-b, --beta-only
- 仅分析和输出 Beta 轨道。
- (注: 这两个选项是互斥的)。

调试 (Debugging)

-v, --verbose
- 启用详细输出模式。
- 将显示操作系统信息、Multiwfn 路径、发送给 Multiwfn 的完整命令序列以及 MOKIT 调用的详细命令。
- 在排查错误时非常有用。

Qchem