QE硅晶体有限q能量损失谱

本篇用硅的两原子原胞，演示有限动量转移下的电子能量损失谱EELS：先求基态，再计算电荷响应，最后提取损失函数和直接/逆介电函数。它与分子吸收光谱使用的扰动、输出和物理含义不同，不能把EELS曲线当作分子吸收强度。

1. 下载与提交顺序#

下载完整结构、公开赝势、真实原生数据和绘图脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd qe-silicon-eels

sbatch job.slurm

CPU分区按56核整节点申请，示例实际使用8个MPI进程；分区维护时等待可用资源。在已经分配的作业内可用 bash run.sh，重新跑整套计算请换一个新的解压目录，避免覆盖已有 work/。

脚本依次执行：一次SCF → 为两个q分别复制同一基态 → 每个q先200步、再续算到600步 → 各自提取有限链与外推链结果 → 独立检查介电函数 → 导出CSV并画图。共13个原生程序分支，两个q使用独立目录，不相互覆盖响应文件。

本例使用鸢算科技集群的QE7.6 CPU构建，Intel oneAPI2024.2编译器、Intel MPI与MKL。绘图及数值提取使用公共 sisl/0.16.4模块提供的Python；用户无需编译或更改公共软件环境。

2. 结构、k点与完整基态输入#

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

&CONTROL calculation='scf',prefix='Si',pseudo_dir='../../pseudo',outdir='./tmp',verbosity='high' /
&SYSTEM ibrav=2,celldm(1)=10.26,nat=2,ntyp=1,ecutwfc=24,ecutrho=96 /
&ELECTRONS conv_thr=1d-10,mixing_beta=.5 /
ATOMIC_SPECIES
Si 28.085 Si.pz-vbc.UPF
ATOMIC_POSITIONS alat
Si 0 0 0
Si .25 .25 .25
K_POINTS automatic
3 3 3 1 1 1

ibrav=2为面心立方原胞，含两个Si原子，celldm(1)=10.26Bohr给出约5.42936Å的常规晶格常数。图中的原胞矢量与原子坐标来自实际QE结构输出，第二个原子按周期等价位置显示。

这里使用QE官方示例随附的NC、LDA-PZ赝势；24/96Ry与3×3×3偏移k网格用于低成本教学流程。改变赝势或采用自己的材料后，需重新检查截断能、k点和电子自洽，不能把本例当作已经完成材料精度收敛的设置。

3. q向量怎么填写#

q是笛卡尔分量，单位为 2π/a，不是倒易晶格基矢的分数坐标，也不是直接填写Å⁻¹。这里的a是前面 celldm(1)对应的晶格参数，不是任一非正交原胞矢量的长度。

例如第一个q的输入：

&LR_INPUT prefix='Si',outdir='./tmp',restart=.false.,restart_step=200 /
&LR_CONTROL itermax=200,q1=0.3,q2=0.,q3=0. /

两个示例分别为 q=(0.3,0,0)和 (0.6,0,0)，实际大小约0.34718与0.69436Å⁻¹。这个流程计算有限q响应，不能直接把q全部设为0来取得光学极限。

当前7.6的输入检查可能把全部分量非正的q误报为 |q| is too small。本例采用正方向；若遇到此提示，不要把非零负向量误认为真实范数接近0。需要其他方向或等价方向对照时，须另外检查实际版本与数值稳定性。本篇不据此给出高精度角向各向异性结论。

续算仍使用同一个q、前置基态和 prefix/outdir：

&LR_INPUT prefix='Si',outdir='./tmp',restart=.true.,restart_step=200 /
&LR_CONTROL itermax=600,q1=0.3,q2=0.,q3=0. /

不要用另一q的断点续算，也不要在两个并发任务里写同一套 tmp/tmp_eels/。

4. 从响应系数得到介电函数#

最终后处理输入为：

&LR_INPUT prefix='Si',outdir='./tmp',eels=.true.,itermax0=600,itermax=10000,extrapolation='osc',units=1,start=0.,end=40.,increment=.04,epsil=0.029399457740524 /

eels=.true.很关键，它按电荷–电荷响应处理EELS。itermax0=600读取600步真实系数，itermax=10000为外推后的后处理维度，不能说实际算了10000步。extrapolation='osc'用于奇偶系数外推；资料包同时保留不外推的200与600步结果供比较。

units=1让能量采样按eV输入；这里0—40eV、间隔0.04eV。epsil仍以Ry输入，0.4eV约为0.0293995Ry，0.8eV约为0.0587989Ry。展宽变大可平滑曲线，但不等于增加了真实能量分辨率或完成了链长收敛。

5. 输出各列与真实结果#

Si.plot_eps.dat的5列依次为：能量、Re(1/ε)、−Im(1/ε)、Re(ε)、Im(ε)。其中 **−Im(1/ε)**是损失函数，是本篇主要绘制的量；不要误取Im(ε)当作EELS。Si.plot_chi.dat另保留复数电荷响应。

峰位对应600步系数、10000维外推、0.4eV展宽；不是实验标定值。本篇用小体系说明操作与q变化，有限链图上的尖峰或平滑程度不能代替k点、截断能和响应步数的收敛比较。

脚本逐能点检查有限链的复数响应，并从原生系数独立组装稀疏矩阵求解；外推链另检查均匀采样点和峰位。直接与逆介电函数逐点互为倒数，电荷响应到逆介电函数的q²、体积和能量换算也逐点检查。这里的数值一致性不等同于材料性质已达到高精度收敛。

6. 导出CSV与重画图片#

先查看资料包中的真实参考结果，不重跑SCF或响应：

bash show-reference.sh
head -n 5 qpos03-eels.csv
head -n 5 qpos06-eels.csv

自己的计算结束后，只重新提取和画图：

bash analyze-results.sh

两个CSV分别对应两个q，表头明确标注损失函数与介电函数。完整复数数组在 eels-fields.npz，报告为 VALIDATION.json，图片在 plots/；原始程序输出在 work/ground/与两个q子目录，参考原件在 reference-results/native/。

参考：官方turboEELS输入说明、光谱后处理输入、QE发行源码硅EELS示例。