Hefei-NAMD：从 H₂ 跃迁偶极到光场驱动

本篇把受驱动两能级模型向前推进一步：先实际计算 H₂ 分子轨道与位置矩阵，再把给定电场转为有效耦合，交给 Hefei-NAMD DEV 做相干传播。

这是固定原子核、冻结两个 RHF 轨道的单电子有效模型。它适合学习“结构→轨道→跃迁矩阵→时变耦合→布居”的操作顺序；结果不是自洽 TDHF/TDDFT 吸收谱，不代表 H₂ 的精确光学激发能，也不是 Hefei-NAMD 的原生电磁场输入功能。

1. 下载，先读取已有结果#

下载输入、量子化学参考、原生传播结果、CSV 和绘图脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd hefei-H2-abinitio-dipole
module load ase/3.29.0
bash read.sh
bash plot.sh
ls analysis-read figures

read.sh 不重新做量子化学或电子传播。两个主要数据表是 H2-abinitio-driven-populations.csv（完整时间序列）和 H2-abinitio-driven-summary.csv（六种条件汇总）。H2-orbital-plane.csv 保存从实际 AO 函数评估得到的轨道截面。

2. 结构与量子化学步骤#

本例结构沿 z 轴，键长 0.740 Å：

H  0  0 -0.370
H  0  0  0.370

采用 PySCF 2.11.0、RHF/STO-3G，SCF 能量收敛 1E-10 Eh。NumPy/SciPy 来自集群 ASE 3.29.0 Python 环境；PySCF 与 h5py 3.14.0 只安装到算例自己的 private-python/。原生传播使用 DEV 2023.03.29 固定快照，ifort 2021.10、Intel MPI 2021.10、MKL 2023.2，选择支持 AVX512 的 CPU 环境。

先检查包内量子化学结果与生成器：

cat reference-results/QUANTUM-MODEL.json
sed -n '1,65p' quantum.py

关键代码如下，int1e_r 的单位为 Bohr：

mol = gto.M(atom='H 0 0 -0.37; H 0 0 0.37',
            unit='Angstrom', basis='sto-3g')
mf = scf.RHF(mol)
mf.conv_tol = 1e-10
mf.conv_tol_grad = 1e-8
mf.kernel()
C = mf.mo_coeff
R = mol.intor_symmetric('int1e_r', comp=3)
position_A = np.einsum('pi,xpq,qj->xij', C, R, C) * 0.529177210903

这里用的是位置矩阵 r01；电子偶极为 −e r01。本例 RHF 总能量 −1.116759307396 Eh，两个轨道能量 −0.578553859832、0.671143491914 Eh，轨道间隔 34.005997329 eV。|z01|=0.492429204 Å，而 x、y 方向为零。正负号随轨道相位而变，不改变同一约定下的布居。

若要亲自重算这一小分子的 SCF，在申请到的 CPU 作业环境、算例根目录运行：

bash setup-quantum.sh

脚本会建立 quantum-new/，重新生成 QUANTUM-MODEL.json、六组 inp/EIGTXT/NATXT/INICON。首次需要下载固定版本的辅助包；不修改公共 Python。不要在登录节点直接运行自己的大分子 SCF。仅复现本篇传播则直接使用下一步，不必重算 SCF。

3. 电场怎样转为 NATXT#

峰值电场 0.1 V/Å，高斯中心 15 fs、宽度 σ=5 fs，演化 40 fs。平行场的光子能量取本例轨道间隔 34.006 eV，属于真空紫外范围。

E(t) = 0.1 exp[-(t-15)^2/(2*5^2)] cos(Eph*t/hbar)  [V/Å]
g(t) = z01 * E(t)                                        [eV]
hbar = 0.6582119281559802                                [eV fs]

单电子长度规范的相互作用是 +e E·r。为了使用 DEV 支持的实反对称耦合格式，写 D12=g/ℏ、D21=−g/ℏ；NATXT 存的是 2*POTIM*D。此约定对应有效哈密顿量的相位变换 S=diag(1,−i)，布居与实非对角 g(t) 的模型相同。不要把数值 eV、V/Å、1/fs 混用，也不要把本例外场耦合解释成核运动产生的非绝热耦合。

cat inputs/parallel/inp
head -3 inputs/parallel/EIGTXT
head -3 inputs/parallel/NATXT
cat inputs/parallel/INICON

平行场输入：

&NAMDPARA
 BMIN=1, BMAX=2, NSW=20004, POTIM=0.002, TEMP=300
 NSAMPLE=1, NTRAJ=1, NAMDTIME=20001, NELM=8, NPAR=1
 ALGO='FSSH', ALGO_INT=0, LHOLE=.FALSE., LSHP=.FALSE.
 LCPTXT=.TRUE., LSPACE=.FALSE., LBINOUT=.FALSE., DEBUGLEVEL='I'
/

LSHP=.FALSE. 关闭表面跳跃；这里只做相干布居传播，不把温度 300 K 当作光子吸收规则。初态 INICON 为 2 1，从第 2 帧、局部能级 1 开始，以保留耦合插值前一帧。

4. 提交传播并读结果#

在算例根目录：

bash submit.sh
squeue -u "$USER"

脚本创建独立 work-new/，通过 cpu-8173m 申请 56 核整节点，程序实际使用 1 MPI、1 线程。不要用 GPU 分区运行这个 AVX512 CPU 可执行文件。

作业结束后：

bash read.sh work-new analysis-new
head -7 analysis-new/H2-abinitio-driven-summary.csv

默认绘图使用 analysis-read/，要绘制自己新运行的结果：

bash read.sh work-new analysis-read
bash plot.sh

每个参考条件只运行一次；如已存在 work-new/，请换一个解压目录再启动，避免覆盖自己的结果。若改了量子化学结构、场强或耦合，需对新输入重新生成独立积分参考，见偏振与步长操作。

5. 得到什么#

本例平行场末态高轨道布居 0.203300；横向电场和关闭电场均为 0。只有给定外场与轨道模型下的布居含义，不能把 0.203300 直接换成实际实验量子产率或吸收截面。

位置积分与 RHF 接口可参阅 PySCF RHF 官方源码文档；程序版本见 PySCF 2.11.0。