ABACUS 实时TDDFT脉冲实操

下载H₂的22个条件、完整输入、原生结果和绘图脚本，解压进入附件根目录，先读取已有数据：

bash analyze.sh

表和数组在reference/analysis，高清图在reference/figures。需要从头重算时，在新的附件目录提交：

sbatch submit.slurm

模板分配1张P100，先顺序运行各条件，再后处理。实时传播使用GPU，线性响应使用同主机的CPU路径；申请了GPU不代表线性响应获得GPU加速。脚本拒绝覆盖已有work，重算请新解压。这是一套含多项控制的演示，整套耗时长于单个小例子。

环境为ABACUS3.10.1模块abacus/3.10.1-p100，GNU13.3、HPC-X2.50/OpenMPI5.0.10rc2、CUDA12.8，1个MPI进程和1线程，双精度；ASE3.29.0环境提供NumPy/SciPy和绘图依赖。不包含Intel Xeon 8173M性能对照。文中的直接计算命令应在已分配的计算资源中执行，不在登录节点运行DFT。

1 先确定是在传播电子#

实时TDDFT给H₂施加短时电场，随时间传播电子态，读取诱导偶极。本例两个H质量均为1.008，键长约0.74Å，放在8Å周期立方盒中，所有原子移动标志均为0。输入虽写calculation md，原子和晶胞实际固定，这不是普通热运动轨迹，也没有进行结构优化。

四个场强0.0025、0.005、0.01、0.02V/Å，两个场方向，三种传播器和三种时间步分别对照；另有辅助网格、盒尺寸和零场轨迹。它们属于一套实时传播工作流，不是十三种独立物性。

2 准备完整输入#

参考目录h2-rt-reference传播4fs；其他大多数场强/传播器条件为2fs。其完整文件为：

INPUT_PARAMETERS
suffix tutorial
pseudo_dir ../../PP_ORB
orbital_dir ../../PP_ORB
ntype 1
basis_type lcao
nbands 8
ecutwfc 40
scf_thr 1e-7
scf_nmax 100
symmetry -1
gamma_only 0
smearing_method fixed
out_chg 0
precision double
cal_force 0
cal_stress 0
calculation md
esolver_type tddft
device gpu
ks_solver elpa
md_type nve
md_nstep 800
md_dt 0.005
md_tfirst 0
td_propagator 0
td_vext 1
td_stype 0
td_tstart 1
td_tend 800
td_vext_dire 3
td_ttype 0
td_gauss_freq 3.66
td_gauss_phase 0
td_gauss_sigma 0.05
td_gauss_t0 40
td_gauss_amp 0.005
out_efield 1
out_dipole 1

ATOMIC_SPECIES
H 1.008 H_ONCV_PBE-1.0.upf
NUMERICAL_ORBITAL
1_H_gga_100Ry_7au_2s1p.orb
LATTICE_CONSTANT
1.8897261258369282
LATTICE_VECTORS
8 0 0
0 8 0
0 0 8
ATOMIC_POSITIONS
Cartesian_angstrom
H
0
2
4 4 3.63 m 0 0 0
4 4 4.37 m 0 0 0

K_POINTS
0
Gamma
1 1 1 0 0 0

• esolver_type tddft选实时电子传播，md_dt 0.005单位fs。
• td_propagator 0/1/2分别为Crank–Nicolson、四阶Taylor和ETRS。本例各有实际结果，仍需针对自己的体系检查时间步稳定性。
• td_stype 0是长度规范；本例分子位于周期盒中央，不能把分子长度规范的偶极流程直接照搬成周期金属的体电导。
• td_vext_dire 3是z方向，沿H–H键；它与静态外场的方向编号约定不同。
• td_gauss_freq 3.66单位fs⁻¹，td_gauss_sigma 0.05单位fs，td_gauss_amp 0.005单位V/Å。
• td_gauss_t0 40是步数，不是0.2fs这个时间本身；40×0.005fs=0.2fs。修改时间步时，应同时改变t0、总步数和场结束步数，保持同一个物理脉冲。
• scf_thr控制每步自洽，不代替传播时间步控制；out_efield和out_dipole分别输出实际场和电子一阶矩。

3 先运行200步的小例子#

h2-rt-dt01用0.01fs走200步，共2fs；脉冲中心设为第20步，仍为0.2fs。在新的附件目录运行：

module purge
module load abacus/3.10.1-p100
python3 -B setup.py
python3 -B run_case.py h2-rt-dt01
head -n 4 work/cases/h2-rt-dt01/OUT.tutorial/SPIN1_DIPOLE

efield_0.dat第一列直接是fs，第二列是V/Å。SPIN1_DIPOLE第一列是步数，后3列是电子密度的一阶矩，单位e·Bohr，不是已经加上离子电荷的总偶极。固定原子时物理诱导偶极为电子一阶矩变化的负值。

python - <<'PY'
import numpy as np
p=np.loadtxt('reference/cases/h2-rt-reference/OUT.tutorial/SPIN1_DIPOLE')
t=p[:,0]*0.005
induced_z=-(p[:,3]-p[0,3])
print('最后时刻/fs',t[-1],'原始诱导偶极最大值/eBohr',abs(induced_z).max())
PY

4 先检查场强和零场，再看谱#

在线性响应范围，诱导偶极除以输入场强后应近似重合；若差异明显，应同时检查非线性、自洽误差、传播步长和文本输出精度。小场强不一定更好：电子一阶矩本身约15e·Bohr，诱导变化很小，原始输出只有有限有效数字，不能凭插值或补零制造更多精度。

附件实际计算了一条4fs零场轨迹。后处理仅对相同方向、盒尺寸、截断能、传播器和时间步的场强条件减去对应零场数值漂移；其余模型保留各自原始诱导变化，比较时使用同样的原始参考，不把不匹配的零场当精确修正。

继续传播时间与频谱提取。输入说明见官方3.10.1参数文档，偶极含义见该版原生输出实现。