«

用于批量生成Q-Chem计算X-ray的程序

wangsy 发布于 阅读：51 密度泛函

📜 qchem_xray.py 详细使用说明

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

1. 脚本概述

qchem_xray.py 是一个用于自动生成 Q-Chem 量子化学计算输入文件的 Python 脚本，专门针对 X 射线光谱 (XAS) 相关的计算。

该脚本的核心功能包括：

2. 核心逻辑与工作流程

A. 配置文件 (input.ini)

脚本的所有行为都由 input.ini 控制。它分为几个关键部分：

B. 依赖关系 (GS -> 激发态)

这是一个至关重要的流程：

  1. 必须首先运行基态 (GS) 计算 (qchem_xray.py GS)。
  2. 必须等待 GS 计算正常完成 (脚本会检查 .out 文件中是否包含 "Have a nice day.")。
  3. GS 计算完成后，才能运行所有其他依赖于 GS 轨道的计算，包括 FCH, FCH_R, LUMO, TD_XAS, TD_XES 和 CACHE。
  4. Z+1 计算不依赖于 GS，可以与 GS 同时运行。

C. 芯轨道自动检测 (缓存机制)

为了运行 FCH, FCH_R, LUMO, TD_XAS, TD_XES 和 CACHE 的计算，脚本必须知道要激发哪个芯轨道 (例如，N 1s)。这个过程是自动的：

  1. 当运行激发态计算（或 CACHE）时，脚本会首先在 [num]/GS/ 目录中寻找 [num]GS_orbital.cache 缓存文件。
  2. 如果缓存文件不存在，或者脚本发现缓存中的轨道局域化程度很差（默认阈值 \< 95.0%），它将自动调用 extract_orbital_composition.py 脚本。
  3. extract_orbital_composition.py 脚本会读取 [num]/GS/ 目录下的 .fchk 文件，生成[num]GS_orbital.cache 文件（未找到 extract_orbital_composition.py 的情况）， 局域化并覆盖当前[num]GS_orbital.cache 文件（找到 extract_orbital_composition.py，但是相应芯轨道阈值 \< 95.0%的情况）。extract_orbital_composition.py 的更多用法请查阅其手册。
  4. 脚本随后读取这个缓存文件，根据 [excited] -> ionization_spin (alpha 或 beta) 的设置，找到能量最低的芯轨道及其编号 (eo)。
  5. 这个轨道编号 (eo) 将被用于生成激发态计算的输入文件 (写入 $occupied 、$reorder_mo或者$alist 部分)。

  6. [num]GS_orbital.cache 文件内容示例如下（此格式由 extract_orbital_composition.py 生成并由 qchem_xray.py 读取）：

    n_alpha=21       n_beta=21
2|98.83              2|98.83            2|49.63          2|49.63
9|54.20              9|54.20            3|49.30          3|49.30
19|45.25            19|45.25

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

    此文件包括局域化后的结果，如果原本波函数的芯轨道就是局域化的，就没有最右两列和最后四行内容，如下所示：

    n_alpha=21       n_beta=21
2|98.83              2|98.83
9|54.20              9|54.20
19|45.25            19|45.25

3. 计算类型与选项详解

A. 支持的计算类型

脚本支持以下 8 种计算类型，必须在命令行中指定一种：

  1. GS (Ground State): 基态计算。这是后续所有激发态计算的基础。
  2. FCH (Final Core-Hole): 芯激发态计算（芯电离）。通过在 $occupied 部分移除一个芯轨道来实现。
  3. FCH_R (FCH Reorganization): 芯激发态重排计算。通过 $reorder_mo 关键字实现，将芯轨道放到最外层轨道，然后控制电荷数来控制最外层电子电离（芯轨道）（通常用于 ROKS/ROHF 计算）。
  4. LUMO (LUMO Excitation): 芯激发到 LUMO 或 LUMO+N 的近似计算。通过在 $occupied 部分将电子从芯轨道移动到指定的未占据轨道。
  5. TD_XAS (Time-Dependent DFT): TDDFT 计算，用于模拟 XAS 吸收谱。通过 $alist 关键字指定芯轨道。
  6. TD_XES (TDDFT for XES): TDDFT 计算，用于模拟 XES 发射谱。通过 $alist 关键字指定芯轨道，内置了FCH的计算。
  7. Z+1 (Z+1 Approximation): Z+1 近似计算。将目标原子替换为周期表中的下一个元素（例如 N -> O）来模拟芯激发。此计算不依赖 GS。
  8. CACHE (Cache Only): 仅生成或更新轨道缓存。用于检查 GS 计算的轨道局域化情况，或在运行激发态计算前强制生成缓存。

B. 命令行选项

4. 完整运行示例

以下是使用示例 input.ini 和 pyridine.xyz 文件运行计算的完整分步指南。

A. 准备文件

  1. 创建一个新的工作目录。
  2. 在该目录中创建以下两个文件：

文件 1: input.ini

# ====================================================================
#  重要提示:
#  此配置文件中的所有值都必须由三引号 """ ... """ 包裹。
#  这是为了确保程序能正确读取，请不要删除或修改这些引号。
#  您只需要修改引号内部的内容即可。
# ====================================================================

# ==========================================================
#               Q-Chem 作业提交配置
# ==========================================================
[submit]
# 设置用于提交计算的命令，如果留空或加#则不自动提交
qsub = """
qqchem 1 2 old
"""

# ==========================================================
#               分子结构与电荷/多重度配置
# ==========================================================
[xyz]
# xyz坐标文件名
filename          = """
pyridine.xyz
"""

# 1. 基态 (GS) 的电荷和自旋多重度
GS_charge_mult    = """
0 1
"""

# 2. 芯电离 (FCH / UKS) 的电荷和自旋多重度
FCH_charge_mult   = """
1 2
"""

# 3. 芯电离 (FCH_R / ROKS) 的电荷和自旋多重度
FCH_R_charge_mult = """
1 2
"""

# 4. LUMO+N 芯激发 (LUMO) 的电荷和自旋多重度 (中性激发)
LUMO_charge_mult  = """
0 1
"""

# 5. 芯激发 (TDDFT/XAS) 的电荷和自旋多重度
TD_XAS_charge_mult   = """
0 1
"""

# 6. XES 发射谱 (TD_XES) 的电荷和自旋多重度 (与FCH相同)
TD_XES_charge_mult = """
1 2
"""

# 7. Z+1 近似的电荷和自旋多重度
Z+1_charge_mult   = """
1 1
"""

# ==========================================================
#               Q-Chem 计算关键词 ($rem 部分)
# ==========================================================
[keys]
# ----------------- 1. 基态 (GS) 计算关键词 -----------------
GS = """
METHOD bp86
BASIS mixed
PURECART 1111
GEN_SCFMAN true
SCF_ALGORITHM diis
BASIS2 def2-svp
UNRESTRICTED true
SCF_GUESS_MIX 5
SCF_CONVERGENCE 8
INTEGRAL_SYMMETRY false
POINT_GROUP_SYMMETRY false
IQMOL_FCHK true
MAX_SCF_CYCLES 300
NBO true
"""

# ----------------- 2. 芯电离 (FCH / UKS) 计算关键词 -----------------
FCH = """
METHOD bp86
BASIS mixed
PURECART 1111
GEN_SCFMAN false
SCF_ALGORITHM diis
SCF_GUESS read
UNRESTRICTED true
MOM_START 1
MOM_METHOD imom
SCF_CONVERGENCE 6
INTEGRAL_SYMMETRY false
POINT_GROUP_SYMMETRY false
IQMOL_FCHK true
MAX_SCF_CYCLES 300
"""

# ----------------- 3. 芯电离 (FCH_R / ROKS) 计算关键词 -----------------
FCH_R = """
METHOD bp86
BASIS mixed
PURECART 1111
GEN_SCFMAN false
SCF_ALGORITHM sgm
SCF_GUESS read
UNRESTRICTED false
SCF_CONVERGENCE 6
INTEGRAL_SYMMETRY false
POINT_GROUP_SYMMETRY false
IQMOL_FCHK true
MAX_SCF_CYCLES 300
"""

# ----------------- 4. 芯激发 (LUMO+N) 计算关键词 -----------------
LUMO = """
METHOD bp86
BASIS mixed
PURECART 1111
GEN_SCFMAN false
SCF_ALGORITHM diis
SCF_GUESS read
UNRESTRICTED true
MOM_START 1
MOM_METHOD imom
SCF_CONVERGENCE 6
INTEGRAL_SYMMETRY false
POINT_GROUP_SYMMETRY false
IQMOL_FCHK true
MAX_SCF_CYCLES 300
"""

# ----------------- 5. 芯激发 (TDDFT/XAS) 计算关键词 -----------------
TD_XAS = """
EXCHANGE bp86
BASIS mixed
PURECART 1111
SKIP_SCFMAN true
SCF_GUESS read
UNRESTRICTED true
CIS_N_ROOTS 300
CIS_TRIPLETS false
TRNSS true
TRTYPE 3
N_SOL 1
FAST_XAS true
! XAS_EDGE 6
! XAS_SCREEN_LEVEL 1
! REL_SHIFT 6
"""

# ----------------- 6. XES 发射 (TDDFT/XES) 计算关键词 -----------------
TD_XES = """
EXCHANGE bp86
BASIS mixed
PURECART 1111
SCF_GUESS read
UNRESTRICTED true
MOM_START 1
MOM_METHOD imom
SCF_CONVERGENCE 6
CIS_N_ROOTS 5
CIS_TRIPLETS false
"""

# ----------------- 7. Z+1 近似计算关键词 -----------------
Z+1 = """
METHOD bp86
BASIS mixed
PURECART 1111
GEN_SCFMAN true
SCF_ALGORITHM diis
BASIS2 def2-svp
UNRESTRICTED true
SCF_GUESS_MIX 5
SCF_CONVERGENCE 8
INTEGRAL_SYMMETRY false
POINT_GROUP_SYMMETRY false
IQMOL_FCHK true
MAX_SCF_CYCLES 300
"""

# ==========================================================
#               基组定义 ($basis 部分)
# ==========================================================
# ----------------- 激发原子配置 -----------------
[excited]
# 定义被激发的原子 (可以用原子序号、范围、元素符号)
# 示例: "1 3-5 N" 表示第1个原子、第3到5个原子、以及所有的N原子
num   = """
N 2
"""

# 为激发原子指定的基组
basis = """
6-311g
"""

# 指定从 'alpha' 还是 'beta' 轨道激发 (默认为 'beta')
ionization_spin = """
beta
"""

# ----------------- 非激发原子配置 -----------------
# 可以定义多组非激发原子及其基组
# 命名规则: num1/basis1, num2/basis2, ...
[non_excited]
num1    = """
1
"""
basis1  = """
3-21G
"""

num2    = """
C
"""
basis2  = """
6-31G
"""

num3    = """
H
"""
basis3  = """
STO-3G
"""

# 如果不需要第4组，可以将其留空
num4    = """
"""
basis4  = """
"""

文件 2: pyridine.xyz

11

N     -1.21390500   -0.76688800    0.00000000
C     -1.23050600    0.58519600    0.00000000
C     -0.07028600    1.38150500    0.00000000
C      1.18362600    0.74775800    0.00000000
C      1.21742000   -0.65680000    0.00000000
C      0.00000000   -1.36257100    0.00000000
H     -2.22120100    1.05147700    0.00000000
H     -0.15352900    2.47192200    0.00000000
H      2.10854100    1.33207700    0.00000000
H      2.16639700   -1.20029000    0.00000000
H     -0.00439300   -2.45750200    0.00000000

B. 步骤 1: 运行基态 (GS) 和 Z+1 计算

[excited] -> num 设置为 "N 2"，脚本会将其解析为 "所有 N 原子 (原子 1)" 和 "原子 2 (C 原子)"。

运行命令:

qchem_xray.py GS
qchem_xray.py Z+1

脚本行为:

  1. 创建目录 1/GS/ 和 2/GS/。
  2. 在 1/GS/ 中生成 1GS.in，在 2/GS/ 中生成 2GS.in。
  3. 创建目录 1/Z+1/ 和 2/Z+1/。
  4. 在 1/Z+1/ 中生成 1Z+1.in (N -> O)，在 2/Z+1/ 中生成 2Z+1.in (C -> N)。
  5. 在每个新创建的目录中执行 qqchem 1 2 old 命令提交作业。

示例: 生成的 1GS.in 文件内容
(注意 BASIS mixed 和混合基组的定义)

$rem
   METHOD bp86
   BASIS mixed
   PURECART 1111
   GEN_SCFMAN true
   SCF_ALGORITHM diis
   BASIS2 def2-svp
   UNRESTRICTED true
   SCF_GUESS_MIX 5
   SCF_CONVERGENCE 8
   INTEGRAL_SYMMETRY false
   POINT_GROUP_SYMMETRY false
   IQMOL_FCHK true
   MAX_SCF_CYCLES 300
   NBO true
$end

$nbo
naomo
$end

$molecule
0 1
 N     -1.21390500   -0.76688800    0.00000000
 C     -1.23050600    0.58519600    0.00000000
 C     -0.07028600    1.38150500    0.00000000
 C      1.18362600    0.74775800    0.00000000
 C      1.21742000   -0.65680000    0.00000000
 C      0.00000000   -1.36257100    0.00000000
 H     -2.22120100    1.05147700    0.00000000
 H     -0.15352900    2.47192200    0.00000000
 H      2.10854100    1.33207700    0.00000000
 H      2.16639700   -1.20029000    0.00000000
 H     -0.00439300   -2.45750200    0.00000000
$end

$basis
N 1
6-311g
****
C 2
6-31G
****
C 3
6-31G
****
C 4
6-31G
****
C 5
6-31G
****
C 6
6-31G
****
H 7
STO-3G
****
H 8
STO-3G
****
H 9
STO-3G
****
H 10
STO-3G
****
H 11
STO-3G
****
$end

注意: non_excited 中 num1=1 (3-21G) 被 num2=C (6-31G) 的规则覆盖了，因为原子 1 (N) 匹配了 excited 规则，而原子 2-6 (C) 匹配了 num2 规则。

注意: 如果计算关键词中含有 NBO，将会自动添加以下模块：

$nbo
naomo
$end

C. 步骤 2: (可选) 检查/生成轨道缓存

注意：如果不运行该步骤，不会生成[num]GS_orbital.cache文件。但是后续运行依赖于GS的计算类型，也会自动调用extract_orbital_composition.py生成[num]GS_orbital.cache文件，因此此步骤不是必选。但是推荐计算完所有GS后，进行此步骤，查看有没有异常轨道或者未局域化轨道，防止造成计算资源浪费。

等待步骤 1 中的 GS 计算全部完成。

运行命令:

qchem_xray.py CACHE

脚本行为:

  1. 脚本开始处理原子 1 (N) 和 2 (C)。
  2. 它将调用 extract_orbital_composition.py (因为缓存尚不存在)。
  3. 该脚本会分析 1/GS/1GS.fchk 和 2/GS/2GS.fchk。
  4. 如果轨道局域化程度低，extract_orbital_composition.py 可能会提示用户是否使用 MOKIT 进行局域化（这是交互式的）。
  5. 最终生成 1/GS/1GS_orbital.cache 和 2/GS/2GS_orbital.cache。
  6. 同时生成一个总结文件 orbital.cache。

D. 步骤 3: 运行激发态计算 ( FCH, FCH_R, LUMO, TD_XAS, TD_XES )

等待 GS 计算完成后，您可以运行所有依赖 GS 的计算。

运行命令:

qchem_xray.py FCH
qchem_xray.py FCH_R
qchem_xray.py LUMO
qchem_xray.py TD_XAS
qchem_xray.py TD_XES

脚本行为 (以 FCH 为例):

  1. 创建 1/FCH/ 和 2/FCH/。
  2. 读取缓存: 脚本读取 1/GS/1GS_orbital.cache 和 2/GS/2GS_orbital.cache (如果缓存不存在，它会像步骤 2 一样自动生成它们)。
  3. 确定轨道: 根据 [excited] -> ionization_spin = beta，它会查找 beta 轨道的最低芯轨道。假设对于原子 1 (N)，它找到的是轨道 2。
  4. 生成输入: 生成 1/FCH/1FCH.in。
  5. 提交作业: 在 1/FCH/ 目录中运行 qqchem 1 2 old。

示例: 生成的 1FCH.in 文件内容
(假设 GS 有 21 个 alpha 和 21 个 beta 电子, 且 N 1s (beta) 是轨道 2)

$occupied
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
1 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
$end

$rem
   METHOD bp86
   BASIS mixed
   PURECART 1111
   GEN_SCFMAN false
   SCF_ALGORITHM diis
   SCF_GUESS read
   UNRESTRICTED true
   MOM_START 1
   MOM_METHOD imom
   SCF_CONVERGENCE 6
   INTEGRAL_SYMMETRY false
   POINT_GROUP_SYMMETRY false
   IQMOL_FCHK true
  *MAX_SCF_CYCLES 300
$end

$molecule
1 2
 N     -1.21390500   -0.76688800    0.00000000
 C     -1.23050600    0.58519600    0.00000000
 C     -0.07028600    1.38150500    0.00000000
 C      1.18362600    0.74775800    0.00000000
 C      1.21742000   -0.65680000    0.00000000
 C      0.00000000   -1.36257100    0.00000000
 H     -2.22120100    1.05147700    0.00000000
 H     -0.15352900    2.47192200    0.00000000
 H      2.10854100    1.33207700    0.00000000
 H      2.16639700   -1.20029000    0.00000000
 H     -0.00439300   -2.45750200    0.00000000
$end

$basis
N 1
6-311g
****
C 2
6-31G
****
C 3
6-31G
****
C 4
6-31G
****
C 5
6-31G
****
C 6
6-31G
****
H 7
STO-3G
****
H 8
STO-3G
****
H 9
STO-3G
****
H 10
STO-3G
****
H 11
STO-3G
****
$end

(注意 $occupied 部分的 beta 轨道 (第二行) 中将会去掉轨道 2, 表示电离的芯轨道是2。)

E. 步骤 4: LUMO 专用额外计算

LUMO 计算的行为与 FCH 类似，但它会将电子移动到未占据轨道，并且受 -l 选项控制。

运行命令 (情况 A: 默认激发到 LUMO):

qchem_xray.py LUMO

运行命令 (情况 B: 激发到 LUMO+2):

qchem_xray.py LUMO -l 2

5. 高级选项详解

A. -r, --retry (重试失败/跳过选项)

示例命令:

# 重试上次未生成的`FCH`计算任务, 读取的是`FCH_skipped.log`文件中的原子列表，可以自己调整。
qchem_xray.py FCH -r"

B. -co ORBITALS, --core-orbital ORBITALS (手动指定芯轨道)

示例命令:
假设 input.ini 中 num = "N 2" (即原子 1 和 2)。自动检测发现原子 1 的芯轨道是 2，但原子 2 (C) 的检测失败了。您手动检查发现原子 2 的芯轨道是 3。您可以使用以下命令：

# 运行 FCH
# - Atom 1: 使用 -co 指定的轨道 2
# - Atom 2: 使用 -co 指定的轨道 3
qchem_xray.py FCH -co "1:2 2:3"

C. 自定义计算类型

Qchem

请先 登录 再评论