QE分子吸收光谱与链长检查

下载水分子完整输入、公开赝势、真实结果、CSV提取与绘图脚本，保存为 example.zip。本例由鸢算科技集群的QE7.6 CPU构建完成：Intel oneAPI2024.2编译器、Intel MPI及MKL；sisl/0.16.4模块提供这里绘图脚本所用的Python环境。

python3 -m zipfile -e example.zip .
cd qe-water-tddfpt

1. 先明确计算什么#

本例计算固定水分子的线性光学响应：pw.x先得到基态电荷密度和波函数，turbo_lanczos.x对x、y、z三个电场方向生成响应系数，turbo_spectrum.x再将系数变成极化率张量和吸收谱。

它不是先逐个求空带再相加的独立粒子吸收，也不是光学带隙、发光谱或激子束缚能的通用替代。这里的 S(E)保持软件原生幅值，张量采用文件表头的单位，不用某个峰的位置或高度直接代表实验水谱。

2. 检查结构、基态输入和赝势#

cat inputs/scf.in
cat PSEUDOPOTENTIAL-SOURCES.json

完整基态输入：

&CONTROL
 calculation='scf', prefix='water', pseudo_dir='../../pseudo',
 outdir='./tmp', tprnfor=.true., verbosity='high'
/
&SYSTEM
 ibrav=1, celldm(1)=16, nat=3, ntyp=2, nbnd=12,
 ecutwfc=60, ecutrho=240
/
&ELECTRONS
 conv_thr=1d-10, mixing_beta=.4, electron_maxstep=150
/
ATOMIC_SPECIES
O 15.999 O_HSCV_PBE-1.0.upf
H 1.008 H_HSCV_PBE-1.0.upf
ATOMIC_POSITIONS angstrom
O 3.175063 3.175063 3.175063
H 4.132263 3.175063 3.175063
H 2.935076 4.102360 3.175063
K_POINTS gamma

celldm(1)=16以Bohr为单位，相当于约8.47Å的立方超胞；原子坐标以Å为单位。ecutwfc/ecutrho=60/240以Ry为单位；conv_thr=1d-10控制电子自洽收敛。本例只做固定几何的响应，没有包含结构优化。水分子有8个价电子、4条双占据能带；nbnd=12还为后面的Davidson初始试探态与跃迁分析提供空带。

两份NC-PBE赝势从同一QE发行源码示例取得，原文件来源和SHA256都在资料包。更换元素、泛函或赝势后，应重新做截断能、真空尺寸和基态收敛检查，不能混用不对应的前置波函数。

3. 一次提交，依次完成计算和提取#

sbatch job.slurm

CPU分区按56核整节点申请，示例主要使用16个MPI进程，单独的Davidson阶段使用4个。分区维护时等待可用资源；在已分配的CPU作业目录内也可运行 bash run.sh。重新计算整套实例请使用新的解压目录，脚本拒绝覆盖已有 work/。

工作流顺序是：SCF → 三个方向各200步Lanczos → 从断点续到400、1600、6400步 → 六个有限链后处理 → 六个外推链/展宽后处理 → 一个指定能区Davidson激发态 → 提取CSV和图片。原生计算输出在 work/water/，完整张量数组在 spectra.npz，最终表格为 absorption.csv。

lanczos-200.in写为：

&LR_INPUT prefix='water', outdir='./tmp', restart=.false., restart_step=100 /
&LR_CONTROL itermax=200, ipol=4 /

续算文件仍使用同一个 prefix/outdir，例如：

&LR_INPUT prefix='water',outdir='./tmp',restart=.true.,restart_step=800 /
&LR_CONTROL itermax=6400,ipol=4 /

ipol=4依次做三个方向；不是四个方向，也不是申请四张卡。restart=.true.读取已有断点，必须保留同一套 tmp/、基态和响应文件。不要把其他结构或泛函的断点拷进来续算。

4. 光谱采样、展宽和外推是三个不同设置#

最终后处理输入为：

&LR_INPUT prefix='water',outdir='./tmp',itermax0=6400,itermax=10000,extrapolation='osc',units=1,
start=0.,end=25.,increment=.02,epsil=0.0110247966526965,ipol=4 /

itermax0=6400表示读取6400步真实计算的系数；itermax=10000表示外推后的后处理矩阵维度，不代表实际做了10000步响应计算。extrapolation='osc'为官方的奇偶系数外推方案，不能靠外推代替原始链长检查。

units=1使光子能量的起点、终点和间隔按eV输入，这里0—25eV、间隔0.02eV。一个容易弄错的地方是：epsil仍以Ry输入；0.15eV对应约0.0110248Ry，0.30eV对应约0.0220496Ry。把 epsil=0.15写进去会变成约2.04eV展宽，谱峰明显被抹平。

5. 实际链长比较与图片#

下表比较相同0.15eV展宽、相同10000维外推方案的整段曲线，与6400步结果的相对L₂差；它不是实验误差。

100、200、400步的未外推谱差别很大，不能作为收敛结果。即使低能峰已经稳定，高能区也可能仍有明显变化；使用哪个能量区间，就应检查那个区间的谱峰、强度和链长变化。这里还未声称截断能、真空尺寸或完整高能光谱都已收敛。

6. 查看参考结果，导出和重新画图#

无需重新计算，先查看资料包中的真实参考结果：

bash show-reference.sh
head -n 5 absorption.csv

若已经完成自己的计算，只重新提取和绘图：

bash analyze-results.sh

absorption.csv依次给出能量、两种展宽下的原生S(E)、x/y/z方向虚部极化率。完整9个复数张量分量在 spectra.npz，原始文本为每个后处理子目录的 water.plot_chi.dat和 water.plot_S.dat。图片写入 plots/，可直接用于对比。

数值脚本从原生系数独立组装稀疏矩阵求解：有限链逐能点检查全部9个复数张量分量，外推链另外检查均匀采样点及峰位，并检查原生S(E)与张量的关系。数值重建一致用于确认提取与单位处理，没有把它当作材料物理收敛证明。实测QE7.6把能量坐标从Ry切换到eV时，原生S(E)这一列的幅值不会随之缩放，因此图和CSV明确标为原生幅值，不将它直接当作已校准的每eV谱密度，也不把有限能区的积分直接解释为价电子数。

继续阅读指定能区的激发态。参考：Lanczos输入说明、Spectrum输入说明。